Example Data

The purpose of this analysis is to demonstrate how UAS-SfM aerial point cloud data processed with cloud2trees can be used to inform treatment planning and design. This analysis will demonstrate how to process the spatial tree list following the methods defined in Churchill et al. (2016) for ICO thinning. The objective is to demonstrate how the ICO method can be implemented to provide a cut-retention list based on target stand average density targets for clump proportions.

UAS data collection missions were completed on the Pike-San Isabel National Forest (PSINF) in October 2025 as part of a demonstration to CFRI. Concurrent field data was collected by CFRI crews which will be used in later data validation efforts.

UAS Collection

  • 130 ha flown in 1 day
  • ~5,000 RGB images captured at 2.5 cm resolution
  • Geo-rectified, wall-to-wall imagery
  • Structure-from-Motion (SfM) point cloud with segmented trees
  • Time commitment:
    • 1.5 field days (7 people; but really 3 people since 4 were primarily observers)
    • 1 day of processing (1 person)
  • Using a ~$30K drone flown at 600 ft
    • Processing time will increase with a smaller drone and a 400 ft flying altitude

Field Collection

  • 26 one‑tenth–acre plots
  • Approximately 1 plot per 10–15 acres, depending on stand density
  • Standard forestry metrics with stem‑mapped overstory trees and saplings
  • Time commitment:
    • ~7 field days
    • 3–4 people

load the standard libraries we use to do work

# bread-and-butter
library(tidyverse) # the tidyverse
library(viridis) # viridis colors
library(harrypotter) # hp colors
library(RColorBrewer) # brewer colors
library(scales) # work with number and plot scales
library(latex2exp)

# visualization
library(mapview) # interactive html maps
library(kableExtra) # tables
library(patchwork) # combine plots
library(ggnewscale) # new scale
library(ggrepel) # repel labels

# spatial analysis
library(terra) # raster
library(sf) # simple features
library(dbscan) # spatial clustering

now, define a function to convert metric to imperial units

calc_imperial_units_fn <- function(df) {
  df %>% 
  # convert to imperial units
    dplyr::mutate(
      dplyr::across(
        .cols = tidyselect::ends_with("_cm")
        , ~ .x * 0.394
        , .names = "{.col}_in"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_m")
        , ~ .x * 3.28
        , .names = "{.col}_ft"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_m2_per_ha")
        , ~ .x * 4.359
        , .names = "{.col}_ftac"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_per_ha") & !tidyselect::ends_with("_m2_per_ha")
        , ~ .x * 0.405
        , .names = "{.col}_ac"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_area_ha")
        , ~ .x * 2.471
        , .names = "{.col}_ac"
      )
      , dplyr::across(
        .cols = tidyselect::ends_with("_m2")
        , ~ .x * 10.764
        , .names = "{.col}_ft2"
      )
    ) %>%
    dplyr::rename_with(
      .fn = function(x){dplyr::case_when(
        stringr::str_ends(x,"_cm_in") ~ stringr::str_replace(x,"_cm_in","_in")
        , stringr::str_ends(x,"_m_ft") ~ stringr::str_replace(x,"_m_ft","_ft")
        , stringr::str_ends(x,"_m2_per_ha_ftac") ~ stringr::str_replace(x,"_m2_per_ha_ftac","_ft2_per_ac")
        , stringr::str_ends(x,"_per_ha_ac") ~ stringr::str_replace(x,"_per_ha_ac","_per_ac")
        , stringr::str_ends(x,"_area_ha_ac") ~ stringr::str_replace(x,"_area_ha_ac","_area_ac")
        , stringr::str_ends(x,"_m2_ft2") ~ stringr::str_replace(x,"_m2_ft2","_ft2")
        , TRUE ~ x
      )}
    )
}

Load Data

load data

# where is the processed data from cloud2trees?
input_dir <- "../data/"
# tree top points
treetops_sf_with_dbh <- sf::st_read(file.path(input_dir, "West_Stand_final_detected_tree_points.gpkg"))
## Reading layer `West_Stand_final_detected_tree_points' from data source 
##   `C:\Data\usfs\psinf_utecreek_ico_ex\data\West_Stand_final_detected_tree_points.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 13964 features and 61 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 465870.1 ymin: 4355809 xmax: 466830.4 ymax: 4356255
## Projected CRS: WGS 84 / UTM zone 13N
# what?
treetops_sf_with_dbh %>% dplyr::glimpse()
## Rows: 13,964
## Columns: 62
## $ treeID                    <chr> "57060_466401.1_4355851.6", "27967_466607.1_…
## $ tree_height_m             <dbl> 24.723, 24.350, 23.548, 22.951, 23.115, 22.8…
## $ crown_area_m2             <dbl> 41.8125, 32.8750, 20.5000, 18.5625, 53.0625,…
## $ fia_est_dbh_cm            <dbl> 34.92759, 34.79764, 34.49266, 34.27861, 34.2…
## $ fia_est_dbh_cm_lower      <dbl> 19.41092, 19.13628, 19.17295, 18.98428, 19.0…
## $ fia_est_dbh_cm_upper      <dbl> 56.16955, 56.11775, 56.01807, 55.60062, 55.2…
## $ dbh_cm                    <dbl> 34.92759, 34.79764, 34.49266, 34.27861, 34.2…
## $ is_training_data          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m                     <dbl> 0.3492759, 0.3479764, 0.3449266, 0.3427861, …
## $ radius_m                  <dbl> 0.1746379, 0.1739882, 0.1724633, 0.1713930, …
## $ basal_area_m2             <dbl> 0.09581358, 0.09510194, 0.09344225, 0.092286…
## $ basal_area_ft2            <dbl> 1.0313374, 1.0236773, 1.0058124, 0.9933675, …
## $ ptcld_extracted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ tree_cbh_m                <dbl> 8.388833, 9.573767, 8.006800, 8.593433, 7.82…
## $ is_training_cbh           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ forest_type_group_code    <chr> "220", "220", "220", "220", "220", "220", "2…
## $ forest_type_group         <chr> "Ponderosa pine group", "Ponderosa pine grou…
## $ hardwood_softwood         <chr> "Softwood", "Softwood", "Softwood", "Softwoo…
## $ comp_trees_per_ha         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ max_crown_diam_height_m   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ is_training_hmd           <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ landfire_stand_id         <dbl> 2144, 1468, 1460, 1583, 1974, 1526, 2145, 21…
## $ crown_dia_m               <dbl> 7.296391, 6.469757, 5.108954, 4.861534, 8.21…
## $ crown_length_m            <dbl> 16.33417, 14.77623, 15.54120, 14.35757, 15.2…
## $ crown_volume_m3           <dbl> 455.3149, 323.8458, 212.3964, 177.6749, 540.…
## $ landfire_tree_kg_per_m3   <dbl> 0.1890150, 0.2128839, 0.2036069, 0.2991796, …
## $ landfire_stand_kg_per_m3  <dbl> 0.08, 0.11, 0.11, 0.17, 0.17, 0.17, 0.08, 0.…
## $ landfire_crown_biomass_kg <dbl> 86.06134, 68.94156, 43.24538, 53.15671, 154.…
## $ FIAnum                    <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ REGION                    <chr> "02", "02", "02", "02", "02", "02", "02", "0…
## $ FORESTNUMB                <chr> "12", "12", "12", "12", "12", "12", "12", "1…
## $ DISTRICTNU                <chr> "11", "11", "11", "11", "11", "11", "11", "1…
## $ spcd                      <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ voleq                     <chr> "200FW2W122", "200FW2W122", "200FW2W122", "2…
## $ TCFV                      <dbl> 31.2, 30.5, 28.9, 27.8, 28.0, 27.5, 27.4, 27…
## $ CFV_GRS                   <dbl> 28.9, 28.9, 26.0, 25.6, 26.0, 25.6, 25.6, 25…
## $ VStump                    <dbl> 0.9291809, 0.9232748, 0.9086702, 0.8982098, …
## $ VTip                      <dbl> 0.6295683, 0.5614567, 0.6773489, 0.6989350, …
## $ BFV_GRS                   <dbl> 150, 150, 130, 130, 130, 130, 130, 130, 120,…
## $ BFV_GRS_INTL              <dbl> 165, 165, 140, 140, 140, 140, 140, 140, 130,…
## $ Poly_Name                 <chr> "West", "West", "West", "West", "West", "Wes…
## $ Desc_                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Poly                      <chr> "Boundary Adjusted", "Boundary Adjusted", "B…
## $ CN                        <chr> "c6916fad-74f1-4326-a788-f3c2eaf388b9", "c69…
## $ Area_MtSq                 <dbl> 69250.85, 69250.85, 69250.85, 69250.85, 6925…
## $ Area_Ac                   <dbl> 71.63097, 71.63097, 71.63097, 71.63097, 71.6…
## $ Area_Ha                   <dbl> 6.925085, 6.925085, 6.925085, 6.925085, 6.92…
## $ Perim_M                   <dbl> 1369.383, 1369.383, 1369.383, 1369.383, 1369…
## $ PerimL_M                  <dbl> 817.6369, 817.6369, 817.6369, 817.6369, 817.…
## $ Desc1                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Shape_Leng                <dbl> 2589.825, 2589.825, 2589.825, 2589.825, 2589…
## $ Shape_Area                <dbl> 289880.3, 289880.3, 289880.3, 289880.3, 2898…
## $ size_class                <chr> "Saw Timber (9-14.9 inch DBH)", "Saw Timber …
## $ dbscan_cluster            <chr> "0", "6", "6", "6", "6", "6", "6", "6", "6",…
## $ clump_id                  <chr> "339", "6", "6", "6", "6", "6", "6", "6", "6…
## $ clump_n_trees             <int> 1, 12235, 12235, 12235, 12235, 12235, 12235,…
## $ clump_class               <chr> "Individual", ">15 trees", ">15 trees", ">15…
## $ largest_in_clump          <chr> "Largest", "Largest", "Second Largest", "", …
## $ geom                      <POINT [m]> POINT (466401.1 4355852), POINT (46660…

generate a stand boundary since we don’t have one at the moment and look at where we are in the US

stand_sf <- 
  treetops_sf_with_dbh %>% 
  sf::st_union() %>% 
  sf::st_concave_hull(ratio = 0.08) %>% 
  sf::st_buffer(2) %>% 
  sf::st_as_sf() %>% 
  dplyr::mutate(dummy=1) %>% 
  sf::st_as_sfc()
# where?
mapview::mapview(
  stand_sf
  , color = "black"
  , lwd = 2
  , alpha.regions = 0
  , label = FALSE
  , legend = FALSE
  , popup = FALSE
  , layer.name = "stand boundary"
)

let’s make a base ggplot2 of the stand

plt_stand <- ggplot2::ggplot() +
  ggplot2::geom_sf(data = stand_sf, color = "black", fill = NA) +
  ggplot2::scale_x_continuous(expand = c(0, 0)) +
  ggplot2::scale_y_continuous(expand = c(0, 0)) +
  ggplot2::labs(
    x = ""
    , y = ""
  ) +
  ggplot2::theme_void() +
  ggplot2::theme(
    legend.position = "top" # c(0.5,1)
    , legend.direction = "horizontal"
    , legend.margin = ggplot2::margin(0,0,0,0)
    , legend.text = ggplot2::element_text(size = 8)
    , legend.title = ggplot2::element_text(size = 8)
    , legend.key = ggplot2::element_rect(fill = NA, color = NA)
    # , plot.title = ggtext::element_markdown(size = 10, hjust = 0.5)
    , plot.title = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
    , plot.subtitle = ggplot2::element_text(size = 8, hjust = 0.5, face = "italic")
  )
# plt_stand

Define Overstory and Clump Spacing

let’s define our DBH limit for classifying trees as “overstory” and we’ll define our distance in meters to separate tree clump groupings

See Churchill et al. (2016) Figure 4 (p.10) and Matonis and Binkley (2018) who “calculated coverage of mosaic-meadows (percentage of stand > 6 m from overstory trees)” (p. 124)

my_clump_dist_m <- 6
my_ostory_dbh_cm <- 5*2.54 # = inches*2.54

plot the stand with some trees

plt_stand +
  ggplot2::geom_sf(
    data = treetops_sf_with_dbh %>% 
      sf::st_intersection(stand_sf) %>% 
      # overstory
      dplyr::mutate(overstory = dbh_cm>=my_ostory_dbh_cm) %>% 
      dplyr::arrange(overstory)
    , mapping = ggplot2::aes(color = overstory)
    , size = 0.8, alpha = 0.7
  ) +
  ggplot2::scale_color_manual(values=c("gray","navy"))

ggplot2::ggsave("../data/01_tree_tops.jpg", width = 7, height = 5)

Define Clumps

create a general function to apply dbscan::dbscan to point data of class sf

With respect to clump size groupings, Churchill et al. (2016) note that:

Proportions for clump sizes should be lumped into four or five bins for operational simplicity. We use 4 or 5 bins (Fig 5): individual trees, small clumps (2-4 trees), medium clumps (5-9 trees), and large clumps (10-20+ trees). Note that when instructed to leave a large clump (e.g. 10-20 trees), marking crews often have difficulty leaving the upper end of the size range (e.g. an 18, 19, or 20 tree clump). Thus adding a fifth bin for “super clumps” may be necessary (e.g. 15-20 trees or 20-25+ trees), especially if the upper size range of clumps is desired. (p. 12-13)

In their prescription worksheet Excel file Churchill et al. (2016) have the following group sizes: 1, 2-4 (Mean clump size: 3), 5-9 (7), 10-15 (12), 16-25 (20). To be able to cut clumps down to achieve the desired largest “16-25” group, we’ll add a “super” group with >25 trees. This same principle can be applied if you desire a clump size group larger than this (i.e. 25-35), whereby a group above this should be added so that the resulting prescription cuts the “super” group to ensure the proper sizing of the desired largest group.

# function to clump data that are sf points
st_clump_points <- function(
  x # point data of class `sf` 
  , clump_dist_m = 6 # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
  , clump_breaks = c(0,1,4,9,15,25,Inf) # where to break the clump size groups
  , clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees") # what to name the clump size groups
) {
  if(!inherits(x,"sf")){stop("`x` requires sf data")}
  if(
    is.na(clump_dist_m) || 
    is.null(clump_dist_m) || 
    !inherits(as.numeric(clump_dist_m),"numeric")
  ){
    stop("`clump_dist_m` must be numeric and not missing")
  }else{clump_dist_m <- as.numeric(clump_dist_m)[1]}
  if(
    any(is.na(clump_breaks)) || 
    any(is.null(clump_breaks)) || 
    !inherits(as.numeric(clump_breaks),"numeric")
  ){
    stop("`clump_breaks` must be numeric and not missing")
  }else{clump_breaks <- as.numeric(clump_breaks)}
  # get points as x,y
    point_clusters <- x %>% 
      dplyr::mutate(
        X = sf::st_coordinates(.)[,1] %>% as.numeric()
        , Y = sf::st_coordinates(.)[,2] %>% as.numeric()
      ) %>% 
      # get rid of columns we'll create
      dplyr::select( -dplyr::any_of(c(
        "hey_xxxxxxxxxx"
        , "dbscan_cluster"
        , "clump_id"
        , "clump_n_trees"
        , "clump_n_trees_grp"
        , "clump_dist_m"
        , "clump_n_trees_lb"
        , "clump_n_trees_ub"
      )))
    #############################################################################
    ##### Identify clusters in each stem map plot                           #####
    #############################################################################
    ### Place trees into clusters using an inter-tree distance of 6 m
    my_dbscan_temp <- point_clusters %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(X,Y) %>% 
      dbscan::dbscan(eps = clump_dist_m, minPts = 2)
    
    # my_dbscan_temp %>% str()
    
    ### append cluster ID to tree points
    point_clusters$dbscan_cluster <- my_dbscan_temp$cluster
    # point_clusters$cluster %>% summary()
    # point_clusters %>% sf::st_drop_geometry() %>% dplyr::count(cluster) %>% dplyr::arrange(desc(n)) %>% dplyr::slice_head(n=11)
  
  # clump_labels
  if( !inherits(clump_labels,"character") ){
    clump_labels <- NULL
  }
    
  ### cluster metrics
  point_clusters <- point_clusters %>% 
    dplyr::group_by(dbscan_cluster) %>% 
    dplyr::mutate(
      # unique dbscan_cluster for individuals
      clump_id = dplyr::case_when(
        # cluster = 0 from dbscan::dbscan() are outliers (e.g. individual trees)
        dbscan_cluster == 0 ~ max(my_dbscan_temp$cluster)+dplyr::row_number()
        , T ~ dbscan_cluster
      ) %>% 
      factor()
    ) %>% 
    dplyr::group_by(clump_id) %>% 
    dplyr::mutate(
      dbscan_cluster = factor(dbscan_cluster)
      , clump_n_trees = dplyr::n()
      , clump_n_trees_grp = base::cut(
          clump_n_trees
          , breaks = clump_breaks
          , labels = NULL
        )
    ) %>%
    dplyr::ungroup() %>% 
    dplyr::mutate(clump_dist_m = clump_dist_m) %>% 
    dplyr::select(-c(dbscan_cluster))
  
  ### clump factor lb and ub
  point_clusters <- point_clusters %>% 
    dplyr::bind_cols(
      # get the cut upper and lower bound
      tidyr::separate_wider_regex(
        point_clusters
        , cols = clump_n_trees_grp
        , patterns = c(
          "[\\(\\[]" 
          , clump_n_trees_lb = ".*"
          , ","
          , clump_n_trees_ub = ".*"
          , "[\\)\\]]"
        )
        , cols_remove = FALSE
      ) %>%
      # "Inf" string becomes numeric Inf here
      dplyr::mutate(
        # dplyr::across(
        #   c(clump_n_trees_lb, clump_n_trees_ub)
        #   , as.numeric
        # )
        # , clump_n_trees_lb = clump_n_trees_lb+1
        # store as factor so we can use cut levels later
        clump_n_trees_lb = factor(
          as.numeric(clump_n_trees_lb)+1
          , levels = as.character( clump_breaks[-length(clump_breaks)]+1 ) ## all but the last
          # , ordered = T
        )
        , clump_n_trees_ub = factor(
          clump_n_trees_ub
          , levels = as.character(clump_breaks[-1]) ## all but the first
          # , ordered = T
        )
      ) %>% 
      dplyr::select(clump_n_trees_lb, clump_n_trees_ub)
    )
  
  ### clump factor
  if(!any(is.null(clump_labels))){
    point_clusters <- point_clusters %>% 
      dplyr::mutate(
        clump_n_trees_grp = clump_n_trees_grp %>% 
          forcats::fct_relabel(
            ~ clump_labels[match(.x, levels(clump_n_trees_grp))]
          )
      )
  }
  
  # ensure final clump factor is ordered ... nah, this messes up when levels aren't in data
  # point_clusters <- 
  #   point_clusters %>% 
  #   dplyr::mutate(
  #     clump_n_trees_grp = clump_n_trees_grp %>% 
  #       # order by the lower bound just to be sure
  #       forcats::fct_reorder(as.numeric(as.character(clump_n_trees_lb))) %>% 
  #       # ordered so can make >,< comparisions
  #       ordered()
  #   )
  
  # return
  return(point_clusters)
}

# st_clump_points(
#   treetops_sf_with_dbh %>% dplyr::slice_sample(n=1111)
#   # treetops_clustered
#   , clump_dist_m = 6
#   , clump_breaks = c(0,1,9,Inf)
#   , clump_labels = NULL
# ) %>%
# # dplyr::glimpse()
# sf::st_drop_geometry() %>%
# dplyr::count(clump_n_trees_grp,clump_n_trees_lb,clump_n_trees_ub)

function to filter the tree points attached to stand bounds and apply the st_clump_points function

# create function to pass a unit id and return list of trees with clump groupings
get_tree_clumps <- function(
    trees_sf
    , stand_sf
    , ostory_ht_m = NA
    , ostory_dbh_cm = NA
    , clump_dist_m = NA
    , clump_breaks = NA
    , clump_labels = NA
){
  if(
    is.na(clump_dist_m) || 
    is.null(clump_dist_m) || 
    !inherits(as.numeric(clump_dist_m),"numeric")
  ){
    stop("`clump_dist_m` must be numeric and not missing")
  }else{clump_dist_m <- as.numeric(clump_dist_m)[1]}
  if(
    any(is.na(clump_breaks)) || 
    any(is.null(clump_breaks)) || 
    !inherits(as.numeric(clump_breaks),"numeric")
  ){
    stop("`clump_breaks` must be numeric and not missing")
  }else{clump_breaks <- as.numeric(clump_breaks)}
  # check data
  if(!inherits(stand_sf,"sf") && !inherits(stand_sf,"sfc")){stop("`stand_sf` requires polygon sf data")}
  if(
    !all( sf::st_is(stand_sf, c("POLYGON","MULTIPOLYGON")) )
  ){
    stop("`stand_sf` requires POLYGON sf data")
  }else{
    stand_sf <- sf::st_union(stand_sf)
  }
  if(!inherits(trees_sf,"sf")){stop("`trees_sf` requires polygon sf data")}
  if(
    !all( sf::st_is(trees_sf, c("POINT","MULTIPOINT")) )
  ){
    stop("`trees_sf` requires POINT sf data")
  }
  
  # check ostory definition
  if(is.na(as.numeric(ostory_dbh_cm)) && is.na(as.numeric(ostory_ht_m))){
    warning("`ostory_dbh_cm` and `ostory_ht_m` are not set...using `ostory_dbh_cm` = 12.7")
    if(
      !any(stringr::str_equal(names(trees_sf), "dbh_cm"))
    ){
      stop("`trees_sf` requires 'dbh_cm' column")
    }
    ostory_dbh_cm <- 5*2.54
    # filter data
    ttops_temp <- trees_sf %>%
      sf::st_intersection(stand_sf %>% sf::st_transform(sf::st_crs(trees_sf))) %>% 
      dplyr::filter(
        dbh_cm>=as.numeric(ostory_dbh_cm)
      )
  }else if(!is.na(as.numeric(ostory_dbh_cm))){
    if(
      !any(stringr::str_equal(names(trees_sf), "dbh_cm"))
    ){
      stop("`trees_sf` requires 'dbh_cm' column")
    }
    # filter data
    ttops_temp <- trees_sf %>% 
      sf::st_intersection(stand_sf %>% sf::st_transform(sf::st_crs(trees_sf))) %>% 
      dplyr::filter(
        dbh_cm>=as.numeric(ostory_dbh_cm)
      )
  }else{
    if(
      !any(stringr::str_equal(names(trees_sf), "tree_height_m"))
    ){
      stop("`trees_sf` requires 'tree_height_m' column")
    }
    # filter data
    ttops_temp <- trees_sf %>% 
      sf::st_intersection(stand_sf %>% sf::st_transform(sf::st_crs(trees_sf))) %>% 
      dplyr::filter(
        tree_height_m>=as.numeric(ostory_ht_m)
      )
  }
  
  ### Place trees into clusters using an inter-tree distance
  ttops_temp <- st_clump_points(
    x = ttops_temp
    , clump_dist_m = clump_dist_m
    , clump_breaks = clump_breaks
    , clump_labels = clump_labels
  )
  
  # get rid of columns we'll create
  ttops_temp <-
    ttops_temp %>% 
    dplyr::select( -dplyr::any_of(c(
      "hey_xxxxxxxxxx"
      , "clump_nn_dist_m"
      , "clump_ostory_ht_m"
      , "clump_ostory_dbh_cm"
    )))
  
  # add distance to nearest tree within clump
  #### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!GW: update to use kdtree:
  #### !!!!!!!!!!!!!!!!!!! the `RANN` [package](https://jefferislab.github.io/RANN/) 
  #### !!!!!!!!!!!!!!!!!!! enables KD-tree processing using X and Y coordinates 
  #### !!!!!!!!!!!!!!!!!!! to build the tree which is a super-fast way to find the `x` number of 
  #### !!!!!!!!!!!!!!!!!!! near neighbors for each point in an input dataset (see `RANN::nn2()`).
  ttops_temp <- 
    ttops_temp %>% 
    dplyr::group_by(clump_id) %>%
    tidyr::nest() %>% 
    dplyr::mutate(
      clump_nn_dist_m = purrr::map(data, function(x){
        # get index of nearest neighbor
        i = sf::st_nearest_feature(x)
        # get dist
        d = sf::st_distance(x, x[i,], by_element=TRUE) %>% as.numeric()
        return(d)
      })
    ) %>% 
    tidyr::unnest(cols = c(data, clump_nn_dist_m)) %>% 
    sf::st_set_geometry("geom") %>% # set it cuz it got lost 
    dplyr::ungroup() %>% 
    dplyr::mutate(
      clump_dist_m = clump_dist_m
      # , clump_id_duplicate = clump_id # can use this even after nesting data by clump_id
      , clump_ostory_ht_m = ifelse(is.na(ostory_ht_m), as.numeric(NA), as.numeric(ostory_ht_m))
      , clump_ostory_dbh_cm = ifelse(is.na(ostory_dbh_cm), as.numeric(NA), as.numeric(ostory_dbh_cm))
    ) %>% 
    dplyr::relocate(
      tidyselect::starts_with("clump_")
      , .after = dplyr::last_col()
    )
  # return
  return(ttops_temp)
}

call that get_tree_clumps() function which uses our st_clump_points() function for overstory trees only

# call it
treetops_clustered <- 
  get_tree_clumps(
    trees_sf = treetops_sf_with_dbh # %>% dplyr::slice_sample(n=2222)
    , stand_sf = stand_sf
    , ostory_dbh_cm = my_ostory_dbh_cm
    , clump_dist_m = my_clump_dist_m
    , clump_breaks = c(0,1,4,9,15,25,Inf) # where to break the clump size groups
    , clump_labels = c("Individual","2-4 trees","5-9 trees","10-15 trees","16-25 trees",">25 trees") # what to name the clump size groups
  )

what?

treetops_clustered %>% dplyr::glimpse()
## Rows: 8,808
## Columns: 68
## $ treeID                    <chr> "57060_466401.1_4355851.6", "27967_466607.1_…
## $ tree_height_m             <dbl> 24.723, 24.350, 21.165, 21.191, 20.848, 20.9…
## $ crown_area_m2             <dbl> 41.8125, 32.8750, 19.0625, 31.6875, 29.4375,…
## $ fia_est_dbh_cm            <dbl> 34.92759, 34.79764, 33.62946, 33.37739, 33.3…
## $ fia_est_dbh_cm_lower      <dbl> 19.41092, 19.13628, 18.30126, 18.45933, 18.4…
## $ fia_est_dbh_cm_upper      <dbl> 56.16955, 56.11775, 54.53841, 54.07352, 54.0…
## $ dbh_cm                    <dbl> 34.92759, 34.79764, 33.62946, 33.37739, 33.3…
## $ is_training_data          <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ dbh_m                     <dbl> 0.3492759, 0.3479764, 0.3362946, 0.3337739, …
## $ radius_m                  <dbl> 0.1746379, 0.1739882, 0.1681473, 0.1668870, …
## $ basal_area_m2             <dbl> 0.09581358, 0.09510194, 0.08882384, 0.087497…
## $ basal_area_ft2            <dbl> 1.0313374, 1.0236773, 0.9560998, 0.9418210, …
## $ ptcld_extracted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ ptcld_predicted_dbh_cm    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ tree_cbh_m                <dbl> 8.388833, 9.573767, 8.747933, 9.090733, 8.29…
## $ is_training_cbh           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA…
## $ forest_type_group_code    <chr> "220", "220", "220", "220", "220", "220", "2…
## $ forest_type_group         <chr> "Ponderosa pine group", "Ponderosa pine grou…
## $ hardwood_softwood         <chr> "Softwood", "Softwood", "Softwood", "Softwoo…
## $ comp_trees_per_ha         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_relative_tree_height <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ comp_dist_to_nearest_m    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ max_crown_diam_height_m   <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ is_training_hmd           <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ landfire_stand_id         <dbl> 2144, 1468, 1411, 1468, 1527, 1468, 1412, 14…
## $ crown_dia_m               <dbl> 7.296391, 6.469757, 4.926574, 6.351833, 6.12…
## $ crown_length_m            <dbl> 16.334166, 14.776234, 12.417068, 12.100267, …
## $ crown_volume_m3           <dbl> 455.31488, 323.84579, 157.80023, 255.61813, …
## $ landfire_tree_kg_per_m3   <dbl> 0.1890150, 0.2128839, 0.1841485, 0.2128839, …
## $ landfire_stand_kg_per_m3  <dbl> 0.08, 0.11, 0.11, 0.11, 0.12, 0.11, 0.12, 0.…
## $ landfire_crown_biomass_kg <dbl> 86.06134, 68.94156, 29.05867, 54.41699, 40.1…
## $ FIAnum                    <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ REGION                    <chr> "02", "02", "02", "02", "02", "02", "02", "0…
## $ FORESTNUMB                <chr> "12", "12", "12", "12", "12", "12", "12", "1…
## $ DISTRICTNU                <chr> "11", "11", "11", "11", "11", "11", "11", "1…
## $ spcd                      <int> 122, 122, 122, 122, 122, 122, 122, 122, 122,…
## $ voleq                     <chr> "200FW2W122", "200FW2W122", "200FW2W122", "2…
## $ TCFV                      <dbl> 31.2, 30.5, 25.0, 24.6, 24.2, 24.2, 24.1, 22…
## $ CFV_GRS                   <dbl> 28.9, 28.9, 21.5, 21.5, 21.5, 21.5, 21.5, 20…
## $ VStump                    <dbl> 0.9291809, 0.9232748, 0.8594608, 0.8452151, …
## $ VTip                      <dbl> 0.6295683, 0.5614567, 0.5501549, 0.5516371, …
## $ BFV_GRS                   <dbl> 150, 150, 110, 110, 110, 110, 110, 110, 110,…
## $ BFV_GRS_INTL              <dbl> 165, 165, 125, 125, 125, 125, 125, 120, 120,…
## $ Poly_Name                 <chr> "West", "West", "West", "West", "West", "Wes…
## $ Desc_                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Poly                      <chr> "Boundary Adjusted", "Boundary Adjusted", "B…
## $ CN                        <chr> "c6916fad-74f1-4326-a788-f3c2eaf388b9", "c69…
## $ Area_MtSq                 <dbl> 69250.85, 69250.85, 69250.85, 69250.85, 6925…
## $ Area_Ac                   <dbl> 71.63097, 71.63097, 71.63097, 71.63097, 71.6…
## $ Area_Ha                   <dbl> 6.925085, 6.925085, 6.925085, 6.925085, 6.92…
## $ Perim_M                   <dbl> 1369.383, 1369.383, 1369.383, 1369.383, 1369…
## $ PerimL_M                  <dbl> 817.6369, 817.6369, 817.6369, 817.6369, 817.…
## $ Desc1                     <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ Shape_Leng                <dbl> 2589.825, 2589.825, 2589.825, 2589.825, 2589…
## $ Shape_Area                <dbl> 289880.3, 289880.3, 289880.3, 289880.3, 2898…
## $ size_class                <chr> "Saw Timber (9-14.9 inch DBH)", "Saw Timber …
## $ X                         <dbl> 466401.1, 466607.1, 466600.6, 466608.1, 4666…
## $ Y                         <dbl> 4355852, 4356234, 4356244, 4356218, 4356208,…
## $ geom                      <POINT [m]> POINT (466401.1 4355852), POINT (46660…
## $ clump_id                  <fct> 471, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ clump_n_trees             <int> 1, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 6…
## $ clump_n_trees_grp         <fct> Individual, >25 trees, >25 trees, >25 trees,…
## $ clump_dist_m              <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,…
## $ clump_n_trees_lb          <fct> 1, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 2…
## $ clump_n_trees_ub          <fct> 1, Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, I…
## $ clump_nn_dist_m           <dbl> NA, 4.031129, 4.756574, 3.905125, 4.596194, …
## $ clump_ostory_ht_m         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ clump_ostory_dbh_cm       <dbl> 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12.7, 12…

let’s see a quick summary of the number of trees in each clump size

treetops_clustered %>% 
  sf::st_drop_geometry() %>% 
  dplyr::count(clump_n_trees_grp) %>% 
  dplyr::mutate(
    pct_tot_trees = scales::percent(n/sum(n),accuracy=0.1)
  ) %>% 
  dplyr::rename(trees = n)
## # A tibble: 6 × 3
##   clump_n_trees_grp trees pct_tot_trees
##   <fct>             <int> <chr>        
## 1 Individual          354 4.0%         
## 2 2-4 trees           731 8.3%         
## 3 5-9 trees           607 6.9%         
## 4 10-15 trees         356 4.0%         
## 5 16-25 trees         414 4.7%         
## 6 >25 trees          6346 72.0%

we can quickly look at the tree groups

plt_stand +
  ggplot2::geom_sf(
    data = treetops_clustered %>% 
      # rename the active geometry column to "geometry"
      sf::st_set_geometry("geometry") %>% 
      dplyr::mutate(geometry = sf::st_buffer(geometry, dist = clump_dist_m/2))
    , mapping = ggplot2::aes(fill = clump_id)
    , color = NA
    , show.legend = F
  ) +
  ggplot2::scale_fill_manual(
    values = c(
      viridis::turbo( n = ceiling(length(unique(treetops_clustered$clump_id))/3) ) 
      , viridis::viridis( n = ceiling(length(unique(treetops_clustered$clump_id))/3) )
      , viridis::cividis( n = ceiling(length(unique(treetops_clustered$clump_id))/3) )
    ) %>% sample() 
  ) 

Clump Polygons and Metrics

Churchill et al. (2016) provide instructions for implementing the clump identification (Plotkin et al. 2002) in ArcGIS:

Use the Buffer tool (in the Proximity toolset within the Analysis toolbox) to create a buffer of distance d/2, one half the inter-tree distance, around each point. This quantity d/2 is meant to approximate the crown radius of a “typical” overstory tree. Set the Dissolve Type option to ALL, which dissolves overlapping buffers, creating a reduced set of spatially non-overlapping polygons stored as a multipart polygon feature…Sanchez Meador et al. (2011) provide some useful examples of how clump attributes can be summarized…The method described here can be modified to use measured or modeled crown radii for each tree in place of d/2 (p.36)

# create function to pass a return from st_clump_points/get_tree_clumps and create clump polygons with summary stats
get_clump_summary <- function(get_tree_clumps_ans){
  # get clump_dist_m
  clump_dist_m <- min(get_tree_clumps_ans$clump_dist_m, na.rm = T)
  # create clump polys and summary
  clump_polys_temp <- 
    get_tree_clumps_ans %>% 
      dplyr::ungroup() %>% 
      # rename the active geometry column to "geometry"
      sf::st_set_geometry("geometry") %>% 
      dplyr::mutate(geometry = sf::st_buffer(geometry, dist = clump_dist_m/2)) %>% 
      dplyr::group_by(clump_id, clump_n_trees_grp, clump_n_trees_lb, clump_n_trees_ub) %>%
      dplyr::summarise(
        # union buffered tree points
        geometry = sf::st_union(geometry)
        # summary metrics
        , n_trees = dplyr::n_distinct(treeID)
        , mean_dbh_cm = mean(dbh_cm, na.rm = T)
        , min_dbh_cm = min(dbh_cm, na.rm = T)
        , max_dbh_cm = max(dbh_cm, na.rm = T)
        , mean_tree_height_m = mean(tree_height_m, na.rm = T)
        , min_tree_height_m = min(tree_height_m, na.rm = T)
        , max_tree_height_m = max(tree_height_m, na.rm = T)
        , loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
        , basal_area_m2 = sum(basal_area_m2, na.rm = T)
        , sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
        , .groups = "drop_last"
      ) %>%
      dplyr::ungroup() %>% 
      sf::st_make_valid() %>% 
      dplyr::mutate(
        clump_area_ha = sf::st_area(geometry) %>% as.numeric() %>% `/`(10000)
        , trees_per_ha = (n_trees/clump_area_ha)
        , basal_area_m2_per_ha = (basal_area_m2/clump_area_ha)
        , pct_stand_basal_area = basal_area_m2/sum(basal_area_m2)
        , pct_stand_n_trees = n_trees/sum(n_trees)
        , qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
      ) %>%
      dplyr::select(-c(sum_dbh_cm_sq)) %>% 
    # convert to imperial units
      calc_imperial_units_fn() %>% 
      dplyr::mutate(clump_dist_m = clump_dist_m)
  # calculate distance between clumps
  clump_polys_temp <- clump_polys_temp %>% 
    dplyr::mutate(
      nearest = sf::st_nearest_feature(clump_polys_temp)
      , distance_nearest_clump_m = sf::st_distance(
          clump_polys_temp
          , clump_polys_temp[nearest,]
          , by_element=TRUE
        ) %>% 
        as.numeric()
    ) %>% 
    dplyr::select(-c(nearest))
  # return
  return(clump_polys_temp)
}

get the clump summary data using get_clump_summary()

# get it
clump_polys <- get_clump_summary(get_tree_clumps_ans = treetops_clustered)
# what?
clump_polys %>% dplyr::glimpse()
## Rows: 824
## Columns: 34
## $ clump_id                 <fct> 471, 472, 473, 474, 475, 476, 477, 478, 479, …
## $ clump_n_trees_grp        <fct> Individual, Individual, Individual, Individua…
## $ clump_n_trees_lb         <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ clump_n_trees_ub         <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ geometry                 <POLYGON [m]> POLYGON ((466404.1 4355852,..., POLYG…
## $ n_trees                  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ mean_dbh_cm              <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ min_dbh_cm               <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ max_dbh_cm               <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ mean_tree_height_m       <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ min_tree_height_m        <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ max_tree_height_m        <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ loreys_height_m          <dbl> 24.723, 23.115, 22.135, 21.374, 21.033, 20.72…
## $ basal_area_m2            <dbl> 0.09581358, 0.09204683, 0.09099974, 0.0890104…
## $ clump_area_ha            <dbl> 0.002826142, 0.002826142, 0.002826142, 0.0028…
## $ trees_per_ha             <dbl> 353.8393, 353.8393, 353.8393, 353.8393, 353.8…
## $ basal_area_m2_per_ha     <dbl> 33.90261, 32.56979, 32.19929, 31.49541, 30.56…
## $ pct_stand_basal_area     <dbl> 0.0001779436, 0.0001709480, 0.0001690034, 0.0…
## $ pct_stand_n_trees        <dbl> 0.0001135332, 0.0001135332, 0.0001135332, 0.0…
## $ qmd_cm                   <dbl> 34.92759, 34.23414, 34.03887, 33.66477, 33.16…
## $ mean_dbh_in              <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ min_dbh_in               <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ max_dbh_in               <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ qmd_in                   <dbl> 13.76147, 13.48825, 13.41132, 13.26392, 13.06…
## $ mean_tree_height_ft      <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ min_tree_height_ft       <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ max_tree_height_ft       <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ loreys_height_ft         <dbl> 81.09144, 75.81720, 72.60280, 70.10672, 68.98…
## $ basal_area_ft2_per_ac    <dbl> 147.7815, 141.9717, 140.3567, 137.2885, 133.2…
## $ trees_per_ac             <dbl> 143.3049, 143.3049, 143.3049, 143.3049, 143.3…
## $ clump_area_ac            <dbl> 0.006983396, 0.006983396, 0.006983396, 0.0069…
## $ basal_area_ft2           <dbl> 1.0313374, 0.9907921, 0.9795212, 0.9581088, 0…
## $ clump_dist_m             <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, …
## $ distance_nearest_clump_m <dbl> 0.006544238, 0.267026450, 0.006544237, 1.6041…

we should have the same number of polygons in the summary data as identified in the get_tree_clumps() function

# do these numbers match
identical(
  # clump 
  nrow(clump_polys)
  # clumps in tree list data
  , treetops_clustered %>% dplyr::distinct(clump_id) %>% nrow()
)
## [1] TRUE

plot the clump group size

plt_stand +
  ggplot2::geom_sf(
    data = clump_polys
    , mapping = ggplot2::aes(fill = clump_n_trees_grp)
    , color = NA
  ) +
  harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
  labs(fill = "") +
  ggplot2::guides(size = "none", fill = ggplot2::guide_legend(override.aes = list(size = 5)))

ggplot2::ggsave("../data/02_tree_clumps.jpg", width = 7, height = 5)

Clump Spacing

See Churchill et al. (2016) Figure 4 (p.10) and Matonis and Binkley (2018) who “calculated coverage of mosaic-meadows (percentage of stand > 6 m from overstory trees)” (p. 124)

Since we already buffered the tree points to approximate the crown radius, we’ll continue to use our \(d/2\) where \(d\) is maximum distance between trees for determining tree clumps and is meant to approximate the crown radius of a “typical” overstory tree

# create function to pass a return from get_clump_summary() and get a distance raster
get_clump_dist_rast <- function(get_clump_summary_ans, stand_sf){
  # get clump_dist_m
  clump_dist_m <- min(get_clump_summary_ans$clump_dist_m, na.rm = T)
  # rasterize the clump polygons and then calculate distance between clumps as raster
  dist_rast <- 
    terra::rasterize(
        x = get_clump_summary_ans %>% terra::vect()
        , y = get_clump_summary_ans %>% 
          terra::vect() %>% 
          terra::rast(res = 0.2)
      ) %>% 
      terra::distance() %>% 
      # crop it to stand extent
      terra::crop(
        stand_sf %>% 
          terra::vect()
        , mask = T
      )
  ######### part 2
  # now create openings vector data
  openings_vect <- 
    dist_rast %>% 
      terra::classify(rcl = c(clump_dist_m/2,Inf), others = NA, include.lowest = T) %>% 
      terra::as.polygons(na.rm = T) %>% 
      sf::st_as_sf() %>% 
      sf::st_cast("POLYGON") %>% 
      dplyr::mutate(layer = dplyr::row_number()) %>% 
      dplyr::mutate(
        openining_area_m2 = sf::st_area(geometry) %>% as.numeric()
        , clump_dist_m = clump_dist_m
      )
  
  # return
  return(list(dist_rast = dist_rast, openings_vect = openings_vect))
}
# get it
get_clump_dist_rast_ans <- get_clump_dist_rast(clump_polys, stand_sf = stand_sf)

what is in it?

names(get_clump_dist_rast_ans)
## [1] "dist_rast"     "openings_vect"

just view the distance raster

terra::plot(get_clump_dist_rast_ans$dist_rast, axes = F, main = "Distance to Nearest Clump (m)")

plot the distance raster and openings vector data we just got with overlaid tree clumps and tree points

plt_fnl_1 <- 
  # ggplot2::ggplot() + 
  plt_stand +
    # distance
    ggplot2::geom_tile(
      data = get_clump_dist_rast_ans$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::as.data.frame(xy = T) %>% dplyr::rename(f=3)
      , mapping = ggplot2::aes(x=x, y=y, fill = f)
    ) +
    harrypotter::scale_fill_hp(
      "mischief"
      , na.value = "transparent"
      , name = "distance to\nnearest tree (m)"
      , limits = c(0,22)
    ) +
    # openings
    ggplot2::geom_sf(
      data = get_clump_dist_rast_ans$openings_vect
      , mapping = ggplot2::aes(color = openining_area_m2), fill = NA
    ) +
    ggplot2::scale_color_gradient(
      low = "gray77", high = "gray11"
      , labels = scales::comma_format(accuracy = 1)
      , name = latex2exp::TeX("opening\narea ($\\m^2$)")
      , limits = c(0,11100)
    ) +
    # clumps
    ggnewscale::new_scale_fill() +
    ggplot2::geom_sf(
      data = clump_polys
      , mapping = ggplot2::aes(fill = clump_n_trees_grp)
      , color = NA
    ) +
    harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
    # tree points
    ggplot2::geom_sf(data = treetops_clustered, color = "gray88", shape = ".") +
    ggplot2::theme_void() +
    ggplot2::theme(
      legend.position = "right"
      , legend.direction = "vertical"
      , plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
      , legend.title= ggplot2::element_text(size=8)
      , legend.text= ggplot2::element_text(size = 7)
    )

plot it

# plot
plt_fnl_1

# save it
ggplot2::ggsave("../data/03_clumps_dist.jpg", width = 7, height = 5)

highlight the openings

plt_open_1 <-
  plt_stand +
      # clumps
      ggplot2::geom_sf(
        data = clump_polys
        , mapping = ggplot2::aes(fill = clump_n_trees_grp)
        , color = NA
      ) +
      harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) + 
      # openings
      ggnewscale::new_scale_fill() +
      ggplot2::geom_sf(
        data = get_clump_dist_rast_ans$openings_vect
        , mapping = ggplot2::aes(fill = openining_area_m2), color = NA
      ) +
      ggplot2::scale_fill_gradient(
        low = "gray77", high = "gray11"
        , labels = scales::comma_format(accuracy = 1)
        , name = latex2exp::TeX("opening\narea ($\\m^2$)")
        , limits = c(0,11100)
      ) +
      # tree points
      ggplot2::geom_sf(data = treetops_clustered, color = "gray88", shape = ".") +
      ggplot2::theme_void() +
      ggplot2::theme(
          legend.position = "right"
          , legend.direction = "vertical"
          , plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
          , legend.title= ggplot2::element_text(size=8)
          , legend.text= ggplot2::element_text(size = 7)
        )

# plot
plt_open_1

ggplot2::ggsave("../data/04_clumps_openings.jpg", width = 7, height = 5)

combine them with patchwork

plt_fnl_1 + (plt_open_1 + theme(legend.position = "none")) 

# &
#   ggplot2::theme(
#     plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
#     , legend.title=ggplot2::element_text(size=8)
#     , legend.text=ggplot2::element_text(size = 7)
#   )
ggplot2::ggsave("../data/05_clumps_combined.jpg", width = 10, height = 5.5)

Clump Metrics

create a function to summarize by the tree clump grouping variable (e.g. “Individual”,“2-4 trees”,“5-9 trees”,“10-15 trees”,…)

# create a function to summarize by number of tree clump grouping
get_clump_n_trees_grp_summary = function(get_tree_clumps_ans, get_clump_summary_ans, stand_sf){
  # get area of harvest unit
  #...will use this area in the area calculations such that...
  #...TPA = trees in a certain group size across the whole stand area
  harvest_area_ha <- stand_sf %>%
    sf::st_union() %>% 
    sf::st_as_sf() %>% 
    dplyr::mutate(harvest_area_ha = sf::st_area(.) %>% as.numeric() %>% `/`(10000)) %>% 
    dplyr::pull(harvest_area_ha) %>% 
    .[1]
  # collapse and calculate silv metrics
  dta <- get_tree_clumps_ans %>% 
    sf::st_drop_geometry() %>% 
    dplyr::mutate(stand_area_ha = harvest_area_ha) %>% 
    dplyr::group_by(stand_area_ha,clump_n_trees_grp, clump_n_trees_lb, clump_n_trees_ub) %>%
    dplyr::summarise(
      # summary metrics
      n_trees = dplyr::n_distinct(treeID)
      , mean_dbh_cm = mean(dbh_cm, na.rm = T)
      , mean_tree_height_m = mean(tree_height_m, na.rm = T)
      , loreys_height_m = sum(basal_area_m2*tree_height_m, na.rm = T) / sum(basal_area_m2, na.rm = T)
      , basal_area_m2 = sum(basal_area_m2, na.rm = T)
      , sum_dbh_cm_sq = sum(dbh_cm^2, na.rm = T)
      , .groups = "drop_last"
    ) %>%
    dplyr::ungroup() %>% 
    # attach clump area
    dplyr::left_join(
      get_clump_summary_ans %>% 
        sf::st_drop_geometry() %>% 
        dplyr::group_by(clump_n_trees_grp) %>%
        dplyr::summarise(
          clump_area_ha = sum(clump_area_ha)
          , stand_n_clumps = dplyr::n()
          , .groups = "drop_last"
        ) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by(clump_n_trees_grp)
    ) %>% 
    dplyr::mutate(
      trees_per_ha = (n_trees/stand_area_ha) # (n_trees/clump_area_ha) ... this was not right
      , basal_area_m2_per_ha = (basal_area_m2/stand_area_ha) # (basal_area_m2/clump_area_ha) ... this was not right
      , qmd_cm = sqrt(sum_dbh_cm_sq/n_trees)
      # stand calcs
      , stand_trees_per_ha = sum(n_trees)/stand_area_ha
      , stand_basal_area_m2 = sum(basal_area_m2)
      , stand_basal_area_m2_per_ha = sum(basal_area_m2)/stand_area_ha
      , pct_stand_basal_area = basal_area_m2/stand_basal_area_m2
      , pct_stand_n_trees = n_trees/sum(n_trees)
      , stand_qmd_cm = sqrt(sum(get_tree_clumps_ans$dbh_cm^2, na.rm = T)/sum(n_trees))
    ) %>%
    dplyr::select(-c(sum_dbh_cm_sq)) %>% 
    # convert to imperial units
    calc_imperial_units_fn()
  # return
  return(dta)
}

call it

# call it
clump_n_trees_grp_summary <- get_clump_n_trees_grp_summary(
  get_tree_clumps_ans = treetops_clustered
  , get_clump_summary_ans = clump_polys
  , stand_sf = stand_sf
) 
# what?
clump_n_trees_grp_summary %>% dplyr::glimpse()
## Rows: 6
## Columns: 33
## $ stand_area_ha               <dbl> 30.4578, 30.4578, 30.4578, 30.4578, 30.457…
## $ clump_n_trees_grp           <fct> Individual, 2-4 trees, 5-9 trees, 10-15 tr…
## $ clump_n_trees_lb            <fct> 1, 2, 5, 10, 16, 26
## $ clump_n_trees_ub            <fct> 1, 4, 9, 15, 25, Inf
## $ n_trees                     <int> 354, 731, 607, 356, 414, 6346
## $ mean_dbh_cm                 <dbl> 27.69943, 27.64553, 26.96353, 26.74255, 26…
## $ mean_tree_height_m          <dbl> 15.36712, 15.19598, 14.60823, 14.55307, 14…
## $ loreys_height_m             <dbl> 16.79629, 16.48573, 15.94476, 16.09756, 15…
## $ basal_area_m2               <dbl> 22.22600, 45.55548, 36.09377, 21.00662, 24…
## $ clump_area_ha               <dbl> 1.0004541, 1.7604321, 1.3383869, 0.7875838…
## $ stand_n_clumps              <int> 354, 276, 98, 30, 20, 46
## $ trees_per_ha                <dbl> 11.62264, 24.00042, 19.92921, 11.68830, 13…
## $ basal_area_m2_per_ha        <dbl> 0.7297307, 1.4956915, 1.1850417, 0.6896957…
## $ qmd_cm                      <dbl> 28.27379, 28.16869, 27.51546, 27.40994, 27…
## $ stand_trees_per_ha          <dbl> 289.187, 289.187, 289.187, 289.187, 289.18…
## $ stand_basal_area_m2         <dbl> 538.4492, 538.4492, 538.4492, 538.4492, 53…
## $ stand_basal_area_m2_per_ha  <dbl> 17.67853, 17.67853, 17.67853, 17.67853, 17…
## $ pct_stand_basal_area        <dbl> 0.04127779, 0.08460496, 0.06703281, 0.0390…
## $ pct_stand_n_trees           <dbl> 0.04019074, 0.08299273, 0.06891462, 0.0404…
## $ stand_qmd_cm                <dbl> 27.89901, 27.89901, 27.89901, 27.89901, 27…
## $ mean_dbh_in                 <dbl> 10.91357, 10.89234, 10.62363, 10.53656, 10…
## $ qmd_in                      <dbl> 11.13987, 11.09846, 10.84109, 10.79952, 10…
## $ stand_qmd_in                <dbl> 10.99221, 10.99221, 10.99221, 10.99221, 10…
## $ mean_tree_height_ft         <dbl> 50.40415, 49.84282, 47.91499, 47.73406, 46…
## $ loreys_height_ft            <dbl> 55.09184, 54.07318, 52.29883, 52.80000, 50…
## $ basal_area_ft2_per_ac       <dbl> 3.180896, 6.519719, 5.165597, 3.006383, 3.…
## $ stand_basal_area_ft2_per_ac <dbl> 77.06072, 77.06072, 77.06072, 77.06072, 77…
## $ trees_per_ac                <dbl> 4.707168, 9.720169, 8.071330, 4.733762, 5.…
## $ stand_trees_per_ac          <dbl> 117.1207, 117.1207, 117.1207, 117.1207, 11…
## $ stand_area_ac               <dbl> 75.26123, 75.26123, 75.26123, 75.26123, 75…
## $ clump_area_ac               <dbl> 2.472122, 4.350028, 3.307154, 1.946119, 2.…
## $ basal_area_ft2              <dbl> 239.2406, 490.3592, 388.5133, 226.1152, 25…
## $ stand_basal_area_ft2        <dbl> 5795.868, 5795.868, 5795.868, 5795.868, 57…

summary table

# table it
clump_n_trees_grp_summary %>% 
  dplyr::select(
    clump_n_trees_grp, n_trees
    , mean_dbh_in
    , qmd_in
    , mean_tree_height_ft
    , loreys_height_ft
    , trees_per_ac
    , basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      .cols = c(pct_stand_basal_area, pct_stand_n_trees) 
      , .fns = ~ scales::percent(.x, accuracy = 1)
    )
  ) %>% 
  kableExtra::kbl(
    digits = 1
    , escape = F
    , caption = paste0("Overstory tree clump summary")
    , col.names = c(
      "", "trees"
      , "mean<br>DBH (in)"
      , "QMD (in)"
      , "mean<br>Ht. (ft)"
      , "Loreys<br>Ht. (ft)"
      , "TPA"
      , "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
      , "Pct<br>stand BA"
      , "Pct<br>stand trees"
      )
  ) %>% 
  kableExtra::kable_styling()
Overstory tree clump summary
trees mean
DBH (in)
QMD (in) mean
Ht. (ft)
Loreys
Ht. (ft)
TPA BA
ft2 ac-1
Pct
stand BA
Pct
stand trees
Individual 354 10.9 11.1 50.4 55.1 4.7 3.2 4% 4%
2-4 trees 731 10.9 11.1 49.8 54.1 9.7 6.5 8% 8%
5-9 trees 607 10.6 10.8 47.9 52.3 8.1 5.2 7% 7%
10-15 trees 356 10.5 10.8 47.7 52.8 4.7 3.0 4% 4%
16-25 trees 414 10.5 10.7 46.9 50.8 5.5 3.5 4% 5%
>25 trees 6346 10.9 11.0 48.6 51.4 84.4 55.7 72% 72%

ICO Implementation

Churchill et al. (2016) describe the full process for implementing the ICO approach in The ICO Approach to Quantifying and Restoring Forest Spatial Pattern: Implementation Guide in which the authors lay out the prescription development process:

  1. Identify skips and other special treatment areas
  2. Consider the need for openings
  3. Determine the stand average density target
  4. Determine the appropriate distance to define clumps
  5. Obtain targets for clump proportions
  6. Select target clump proportions for your stand
  7. Generate clump targets for the whole unit
  8. Combine clump and opening targets with leave tree criteria into marking guidelines

The objective here is to: 1) provide the manager with the current conditions (completed above); 2) take the “targets” as set by the manager (steps 3, 5, 6, 7); 3) create the prescription with the leave tree marking.

Let’s implement this prescription development process with our UAS tree list

3. Determine the stand average density target

Step 3 in Churchill et al. (2016):

An average BA, TPA, or SDI target for the stand should be selected that is appropriate for the species, structure, site conditions, and management objectives. Expected mortality from prescribed fire should be factored in. Stand average targets can come from historical reference stands, plant association based stocking guides, density management tools, or a combination of both (see Franklin et al. (2013) for a full discussion of setting density targets). In dry forests, the number and size of old trees must be accounted in setting the density target. To use the ICO method, the target must be converted to TPA (see Table 1). A lower diameter cutoff also needs to be specified for the TPA target. This should be the lower limit in the contract or cutting guidelines given to the marking crew or contractor. (p.11)

this is what Table 1 looks like with TPA values are derived from the formula:

\[ TPA = \frac{BA}{QMD^{2} \times 0.005454} \]

# function to get tpa from ba and qmd
get_tpa <- function(ba_ft2_ac, qmd_in){
  tpa <- round(ba_ft2_ac/((qmd_in^2)*0.005454))
  return(tpa)
} 
# table it
tidyr::crossing(
    ba = seq(40,200,20)
    , qmd = seq(8,20,2)
  ) %>% 
  dplyr::mutate(
    tpa = get_tpa(ba,qmd)
  ) %>% 
  tidyr::pivot_wider(names_from = ba, values_from = tpa) %>% 
  dplyr::mutate(l = "QMD (in)") %>% 
  dplyr::relocate(l) %>% 
  kableExtra::kbl(
    col.names = c(".","", seq(40,200,20))
    , escape = F
    , caption = "Basal Area and QMD to TPA conversion chart"
  ) %>% 
  kableExtra::add_header_above(
    c("","", "Basal Area (ft<sup>2</sup> ac<sup>-1</sup>)"=length(seq(40,200,20)))
    , escape = F
  ) %>%
  kableExtra::kable_styling() %>% 
  kableExtra::column_spec(1:2, bold = T) %>%
  kableExtra::collapse_rows(columns = 1, valign = "middle")
Basal Area and QMD to TPA conversion chart
Basal Area (ft2 ac-1)
. 40 60 80 100 120 140 160 180 200
QMD (in) 8 115 172 229 286 344 401 458 516 573
10 73 110 147 183 220 257 293 330 367
12 51 76 102 127 153 178 204 229 255
14 37 56 75 94 112 131 150 168 187
16 29 43 57 72 86 100 115 129 143
18 23 34 45 57 68 79 91 102 113
20 18 28 37 46 55 64 73 83 92

Current Stand Conditions

For determining targets, the silviculturist needs to know the current conditions. Provide the current stand conditions based on the UAS tree list for the selected stand that are required to set the targets:

  • Current BA
  • Current QMD
  • Current proportion of trees by clump size
clump_n_trees_grp_summary %>% 
  dplyr::select(clump_n_trees_grp, pct_stand_n_trees) %>% 
  dplyr::mutate(
    pct_stand_n_trees = scales::percent(pct_stand_n_trees,accuracy = 1)
  ) %>% 
  tidyr::pivot_wider(names_from = clump_n_trees_grp, values_from = pct_stand_n_trees) %>% 
  kableExtra::kable(
    caption = paste0(
      "Current stand BA (ft<sup>2</sup> ac<sup>-1</sup>): "
      , clump_n_trees_grp_summary$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand QMD (in): "
      , clump_n_trees_grp_summary$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand TPA: "
      , clump_n_trees_grp_summary$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
    )
    , escape = F
    , digits = 1
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::footnote(general = "values are the percent of trees in each clump size")
Current stand BA (ft2 ac-1): 77.1
Current stand QMD (in): 11.0
Current stand TPA: 117
Individual 2-4 trees 5-9 trees 10-15 trees 16-25 trees >25 trees
4% 8% 7% 4% 5% 72%
Note:
values are the percent of trees in each clump size

5. Obtain targets for clump proportions

Step 5 in Churchill et al. (2016):

ICO prescriptions are based on a target proportion of trees in different sized clumps within a stand. Proportions are just the percentage of trees, or TPA, that are in different sized clumps. Basal area proportions can be used, but we have found TPA targets to be more straightforward to use. Ideally, a table summarizing clump proportions for a range of reference conditions in your area is available (Table 2). If not, instructions for developing one are provided in section VI. (p.12)

Section VI of Churchill et al. (2016) notes that

reference spatial information may already be available and summarized in a way that it can be directly incorporated into ICO prescriptions. Such data exist and have been published for areas in Arizona (Abella and Denton 2009, Sánchez Meador et al. 2011), the eastern Washington Cascades (Churchill et al. 2013), the northern Rockies (Larson et al. 2012), and the Sierra Nevada (Lydersen et al. 2013). Reference datasets for using ICO in other forest types, such as coastal Douglas-fir or Pacific silver fir, also exist (Larson and Churchill 2008). (p.28)

Table 2 is:

6. Select target clump proportions for your stand

Now set the desired BA, QMD, and proportion of trees in each clump size:

########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################
# desired BA
  # target_ba_ft2_per_ac <- 9 # cannot be > current BA
  target_ba_ft2_per_ac <- 50 # cannot be > current BA
# desired QMD
  target_qmd_in <- 12
# desired proportion (%) of trees in each clump size
  # !cannot be create larger proportion of ">25 trees" clump as this would require adding trees...
  # c("Individual", "2-4 trees",    "5-9 trees",    "10-15 trees","16-25 trees",">25 trees")
  # c(.18, .33, .24, .10, .15)
  # target_pcts <- c(66, 33, 1, 0, 0, 0)
  target_pcts <- c(25, 40, 25, 10, 0, 0)
########################################################################################
########################################################################################
# desired BA, QMD, and proportion of trees in each clump size
########################################################################################
########################################################################################

Check set up and define data with targets based on Churchill et al. (2016) and their prescription worksheet Excel file to help develop ICO prescriptions

get_target_check_prescription = function(
  clump_n_trees_grp_summary_dta
  , target_ba_ft2_per_ac = NA
  , target_qmd_in = NA
  , target_pcts = NA
){
  if(
    is.na(target_ba_ft2_per_ac) | is.na(target_qmd_in) | any(is.na(target_pcts))
  ){
    stop("must set all of the function parameters:\n`target_ba_ft2_per_ac`, `target_qmd_in`, and `target_pcts`")
  }
  #############################################
  # check target BA and TPA
  #############################################
  if(as.numeric(target_ba_ft2_per_ac)>clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1]){
    stop(
      "target BA in `target_ba_ft2_per_ac` of "
      , round(as.numeric(target_ba_ft2_per_ac),1), " is greater than current BA of "
      , clump_n_trees_grp_summary_dta$stand_basal_area_ft2_per_ac[1] %>% round(1)
    )
  }
  if(
    get_tpa(target_ba_ft2_per_ac, target_qmd_in)>clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
  ){
    stop(
      "target TPA based on `target_ba_ft2_per_ac` and `target_qmd_in` of "
      , round(as.numeric(get_tpa(target_ba_ft2_per_ac, target_qmd_in)),1), " is greater than current TPA of "
      , clump_n_trees_grp_summary_dta$stand_trees_per_ac[1] %>% round(1)
      , "\n adjust `target_ba_ft2_per_ac` and/or `target_qmd_in` to get valid TPA"
    )
  }
  # if(sum(target_pcts)!=1){stop("`target_pcts` must sum to 1 (100%)")}
  #############################################
  # define data with current and target
  # ... this is "smart" in that percentages are adj based on:
  # ... 0) are there missing targets?
  # ... ... if < 6 numbers provided in `target_pcts` then the largest tree groups get targets of 0
  # ... 1) do targets sum to 1? 
  # ... ... if not trees are distributed proportionally based on targets provided and trees available
  # ... 2) is target in largest clump size > current conditions?
  # ... ... if yes, target is set to current condition
  # ... 3) is target listed in clump size > current largest clump with trees?
  # ... ... if yes, target for largest clump size is shifted to current largest clump with trees
  #############################################
    # define clump group data even if levels are missing from data
    clump_n_trees_grp_df <- 
      dplyr::tibble(
        # get all possible labels from the factor levels
        clump_n_trees_grp = levels(clump_n_trees_grp_summary_dta$clump_n_trees_grp)
        , clump_n_trees_lb = levels(clump_n_trees_grp_summary_dta$clump_n_trees_lb)
        , clump_n_trees_ub = levels(clump_n_trees_grp_summary_dta$clump_n_trees_ub)
      ) %>% 
      # ensure they are still factors
      dplyr::mutate(
        dplyr::across(
          dplyr::everything()
          , ~factor(.x, levels = .x)
        )
        # get numeric values
        , min_clump_n_trees = clump_n_trees_lb %>% as.character() %>% as.numeric()
        , max_clump_n_trees = clump_n_trees_ub %>% as.character() %>% as.numeric()
      ) %>% 
      # calculate the median to use as the "average" of the clump size
      dplyr::rowwise() %>% 
      dplyr::mutate(
        mean_clump_n_trees = median(c( # note the c()
          min_clump_n_trees
          , ifelse( is.infinite(max_clump_n_trees), ceiling(min_clump_n_trees*1.3), max_clump_n_trees )
        )) %>% 
        floor()
          # c(1,3,7,12,20,30)
      ) %>% 
      dplyr::ungroup()
    
    target_data <-  
      # create data for joining if missing clump groups
      clump_n_trees_grp_df %>% 
      dplyr::mutate(
        stand_area_ac = clump_n_trees_grp_summary_dta$stand_area_ac[1]
      ) %>% 
      dplyr::left_join(
        clump_n_trees_grp_summary_dta %>% 
          dplyr::ungroup() %>% 
          dplyr::select(clump_n_trees_grp, pct_stand_n_trees, stand_n_clumps)
        , by = "clump_n_trees_grp"
      ) %>% 
      dplyr::mutate(
        pct_stand_n_trees = dplyr::coalesce(pct_stand_n_trees,0)
        , stand_n_clumps = dplyr::coalesce(stand_n_clumps,0)
      ) %>% 
      # add targets
      dplyr::bind_cols(
        pct_stand_n_trees_target = c(
            as.numeric(target_pcts)
            , rep(0,times=nrow(clump_n_trees_grp_df))
          )[1:nrow(clump_n_trees_grp_df)] # pad target with 0's
      ) %>% 
      # adjust target based on difference from 1
      dplyr::mutate(
        pct_stand_n_trees_target = pct_stand_n_trees_target*(1/sum(pct_stand_n_trees_target))
        # largest clump size with trees
        , largest_w_trees = max(
          ifelse(
            dplyr::coalesce(pct_stand_n_trees)>0
            ,min_clump_n_trees
            ,NA
          )
          ,na.rm = T
        )
        , largest_w_trees_target = max(
          ifelse(
            dplyr::coalesce(pct_stand_n_trees_target)>0
            ,min_clump_n_trees
            ,NA
          )
          ,na.rm = T
        )
      ) %>% 
      # move target for largest clump size to the largest current clump size 
      dplyr::mutate(
        pct_stand_n_trees_target = dplyr::case_when(
          min_clump_n_trees==largest_w_trees &
            largest_w_trees_target>largest_w_trees ~ max(
              ifelse(min_clump_n_trees==largest_w_trees_target,pct_stand_n_trees_target,0)
            )
          , T ~ pct_stand_n_trees_target
        )
      ) %>% 
      # adjust target based on current conditions
      dplyr::mutate(
        pct_stand_n_trees_target = dplyr::case_when(
          min_clump_n_trees>largest_w_trees &
            pct_stand_n_trees_target > 0 ~ 0
          , min_clump_n_trees==largest_w_trees &
            pct_stand_n_trees_target > pct_stand_n_trees ~ pct_stand_n_trees
          , T ~ pct_stand_n_trees_target
        )
      ) %>% 
      # finally, re-scale again based on adjustments
      dplyr::mutate(
        pct_stand_n_trees_target = dplyr::case_when(
          min_clump_n_trees==largest_w_trees ~ pct_stand_n_trees_target
          , T ~ pct_stand_n_trees_target * (
            # pct remaining to scale to
            (1-max(ifelse(min_clump_n_trees==largest_w_trees,pct_stand_n_trees_target,0))) /
            # current pct remaining total allocated
            sum(
              ifelse(min_clump_n_trees!=largest_w_trees,pct_stand_n_trees_target,0))
            )
        )
      ) %>% 
      # add other targets
      dplyr::rename(pct_stand_n_trees_current = pct_stand_n_trees) %>% 
      dplyr::mutate(
        stand_trees_per_ac_current = clump_n_trees_grp_summary_dta$stand_trees_per_ac[1]
        , stand_trees_per_ac_target = dplyr::coalesce(get_tpa(target_ba_ft2_per_ac, target_qmd_in),0)
        , trees_per_acre_current = stand_trees_per_ac_current*pct_stand_n_trees_current
        , trees_per_acre_target = stand_trees_per_ac_target*pct_stand_n_trees_target
        , clumps_per_acre_current = trees_per_acre_current/mean_clump_n_trees
        , clumps_per_acre_target = trees_per_acre_target/mean_clump_n_trees
        , stand_n_clumps_current = stand_n_clumps
        , stand_n_clumps_target = (clumps_per_acre_target*stand_area_ac) %>% round(0)
      ) %>% 
      dplyr::select(-c(tidyselect::starts_with("largest_w_trees"), stand_n_clumps))
  # ????  
  # target_data %>% glimpse()
  
    # issue warning about targets
    if(
      !all(
        target_data$pct_stand_n_trees_target == c(
          as.numeric(target_pcts)
          , rep(0,times=nrow(clump_n_trees_grp_df))
        )[1:nrow(clump_n_trees_grp_df)]
      )
    ){
      warning(
        "proportion of trees in each clump size target `target_pcts` adjusted!!!"
        , "\nfrom : ", paste(round(target_pcts,2),collapse = ",")
        , "\nto : ", paste(round(target_data$pct_stand_n_trees_target,2),collapse = ",")
      )
    }
  # return
  return(target_data)
}
# call it
prescription_target_data <- get_target_check_prescription(
  clump_n_trees_grp_summary
  , target_ba_ft2_per_ac = target_ba_ft2_per_ac
  , target_qmd_in = target_qmd_in
  , target_pcts = target_pcts
)
## Warning in get_target_check_prescription(clump_n_trees_grp_summary, target_ba_ft2_per_ac = target_ba_ft2_per_ac, : proportion of trees in each clump size target `target_pcts` adjusted!!!
## from : 25,40,25,10,0,0
## to : 0.25,0.4,0.25,0.1,0,0
# what?
prescription_target_data %>% dplyr::glimpse()
## Rows: 6
## Columns: 17
## $ clump_n_trees_grp          <fct> Individual, 2-4 trees, 5-9 trees, 10-15 tre…
## $ clump_n_trees_lb           <fct> 1, 2, 5, 10, 16, 26
## $ clump_n_trees_ub           <fct> 1, 4, 9, 15, 25, Inf
## $ min_clump_n_trees          <dbl> 1, 2, 5, 10, 16, 26
## $ max_clump_n_trees          <dbl> 1, 4, 9, 15, 25, Inf
## $ mean_clump_n_trees         <dbl> 1, 3, 7, 12, 20, 30
## $ stand_area_ac              <dbl> 75.26123, 75.26123, 75.26123, 75.26123, 75.…
## $ pct_stand_n_trees_current  <dbl> 0.04019074, 0.08299273, 0.06891462, 0.04041…
## $ pct_stand_n_trees_target   <dbl> 0.25, 0.40, 0.25, 0.10, 0.00, 0.00
## $ stand_trees_per_ac_current <dbl> 117.1207, 117.1207, 117.1207, 117.1207, 117…
## $ stand_trees_per_ac_target  <dbl> 64, 64, 64, 64, 64, 64
## $ trees_per_acre_current     <dbl> 4.707168, 9.720169, 8.071330, 4.733762, 5.5…
## $ trees_per_acre_target      <dbl> 16.0, 25.6, 16.0, 6.4, 0.0, 0.0
## $ clumps_per_acre_current    <dbl> 4.7071680, 3.2400563, 1.1530472, 0.3944802,…
## $ clumps_per_acre_target     <dbl> 16.0000000, 8.5333333, 2.2857143, 0.5333333…
## $ stand_n_clumps_current     <dbl> 354, 276, 98, 30, 20, 46
## $ stand_n_clumps_target      <dbl> 1204, 642, 172, 40, 0, 0
# !!!!!!!!!!!!!!!!!!!! need to ensure sufficient trees in larger clump groups if increasing target for lower groups 
### e.g. see the third record:
# $ pct_stand_n_trees_current  <dbl> 0.28, 0.43, 0.19, 0.07, 0.00, 0.01
# $ pct_stand_n_trees_target   <dbl> 0.40, 0.30, 0.28, 0.02, 0.00, 0.00

current vs target

prescription_target_data %>% 
  dplyr::select(
    clump_n_trees_grp
    , tidyselect::starts_with("pct_stand_n_trees")
    , tidyselect::starts_with("trees_per_acre_")
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_stand_n_trees")
      , ~ scales::percent(.x,accuracy = 1)
    )
  ) %>% 
  kableExtra::kable(
    caption = paste0(
      "Current stand (ft<sup>2</sup> ac<sup>-1</sup>): "
      , clump_n_trees_grp_summary$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand QMD (in): "
      , clump_n_trees_grp_summary$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand TPA: "
      , clump_n_trees_grp_summary$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
    )
    , escape = F
    , digits = 1
    , col.names = c(
      "", rep(c("current","target"),2)
    )
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(
    c(" "=1,"% Trees"=2, "TPA"=2)
  )
Current stand (ft2 ac-1): 77.1
Current stand QMD (in): 11.0
Current stand TPA: 117
% Trees
TPA
current target current target
Individual 4% 25% 4.7 16.0
2-4 trees 8% 40% 9.7 25.6
5-9 trees 7% 25% 8.1 16.0
10-15 trees 4% 10% 4.7 6.4
16-25 trees 5% 0% 5.5 0.0
>25 trees 72% 0% 84.4 0.0

8. Combine clump and opening targets with leave tree criteria into marking guidelines

Use our UAS tree list to generate the prescription:

  • start with the largest clump size currently with trees
  • cut trees to the next largest clump size until desired # clumps is reached
  • repeat with each successive clump size through to individual tree selection
  • if possible, cut in same clump until desired proportions are reached to minimize machine time

Function to cut clumps

put the entire process outlined immediately above into a function to use with the tree list to cut each clump and then we’ll select the desired proportion of clumps

define intermediate functions

define intermediate functions

##############################################
# working with sf LINESTRINGS
##############################################
  # first two functions borrowed from https://github.com/metafor-ulaval/ALSroads/blob/main/R/line_tools.R
  ########
  # Get heading of both ends of a line
  ########
    st_ends_heading <- function(line){
      M <- sf::st_coordinates(line)
      i <- c(2, nrow(M) - 1)
      j <- c(1, -1)
      
      headings <- mapply(i, j, FUN = function(i, j) {
        Ax = M[i-j,1]
        Ay = M[i-j,2]
        Bx = M[i,1]
        By = M[i,2]
        atan2(Ay-By, Ax-Bx)*180/pi
      })
      names(headings) <- c("head", "tail")
      return(headings)
    }
  ########
  # extend the line on both ends
  ########
    st_extend_line <- function(line, distance, end = "BOTH"){
      if (!(end %in% c("BOTH", "HEAD", "TAIL")) | length(end) != 1) stop("'end' must be 'BOTH', 'HEAD' or 'TAIL'")
    
      M <- sf::st_coordinates(line)[,-3]
      keep <- !(end == c("TAIL", "HEAD"))
      
      ends <- c(1, nrow(M))[keep]
      headings <- st_ends_heading(line)[keep] / 180 * pi
      distances <- if (length(distance) == 1) rep(distance, 2) else distance[1:2]
      
      M[ends, 1:2] <- M[ends, 1:2] + distances[keep] * c(cos(headings), sin(headings))
      newline <- sf::st_linestring(M)
    
      # If input is sfc_LINESTRING and not sfg_LINESTRING
      if (is.list(line)) newline <- sf::st_sfc(newline, crs = sf::st_crs(line))
      
      return(newline)
    }
  ########
  # pass an sf dataframe of points and return a line between the farthest points
  ########
    st_points_to_line <- function(pts, line_ext=0) {
      # find farthest distance between points
      dist_temp <- sf::st_distance(pts)
      
      # get the points
      f_pts_temp <-  
        pts %>% 
          dplyr::ungroup() %>% 
          dplyr::slice(
            # get the farthest points from distance matrix
            which(dist_temp == max(dist_temp), arr.ind = TRUE)[1,]
          )
      
      # draw a line between the farthest two points
      f_line_temp <- f_pts_temp %>% 
          # convert to linestring
          dplyr::ungroup() %>% 
          dplyr::summarise(n=dplyr::n()) %>%
          sf::st_cast("LINESTRING") %>%
          dplyr::ungroup() %>% 
          dplyr::select(-c(n))
      
      # and apply the line extension
        farthest_line <- st_extend_line(f_line_temp, distance = line_ext)
      # return
        return(farthest_line)
    }
    # st_points_to_line(ttops_temp %>% dplyr::slice_head(prop = .1), line_ext = 6) %>% 
    #   ggplot() + geom_sf() + geom_sf(data = ttops_temp %>% dplyr::slice_head(prop = .1)) + theme_void()

  ########
  # Function to calculate Euclidean distance between 2 points
  ########
    st_euclidean_distance <- function(p1,p2) {
        return(sqrt((p2[1] - p1[1])^2 + (p2[2] - p1[2])^2))
    }

  ########
  # return a line perpendicular to current line
  ########
  ### https://stackoverflow.com/questions/56771058/perpendicular-lines-at-regular-intervals-along-lines-with-multiple-nodes
    # Function to calculate 2 points on a line perpendicular to another defined by 2 points p1,p2
    # For point at interval, which can be a proportion of the segment length, or a constant
    st_perp_line <- function(interval=0.5, my_line, proportion=TRUE) {
      # get end points of line
      p1 <- my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[1,]
      p2 <- my_line %>% sf::st_cast("POINT") %>% sf::st_coordinates() %>% .[2,]
      # get length of line to return equal length line
      line_len <- sf::st_length(my_line) %>% as.numeric() %>% `/`(2)      
      # get crs of line
      my_crs <- sf::st_crs(my_line)
      
      # Calculate x and y distances
      x_len <- p2[1] - p1[1]
      y_len <- p2[2] - p1[2]
      
      # If proportion calculate reference point from tot_length
      if (proportion) {
          point <- c(p1[1]+x_len*interval,p1[2]+y_len*interval)
      }
      # Else use the constant value
      else {
          tot_len <- st_euclidean_distance(p1,p2)
          point <- c(p1[1]+x_len/tot_len*interval,p1[2]+y_len/tot_len*interval)
      }    
      
      # Calculate the x and y distances from reference point to point on line line_len distance away    
      ref_len <- st_euclidean_distance(point,p2)
      xn_len <- (line_len / ref_len) * (p2[1] - point[1])
      yn_len <- (line_len / ref_len) * (p2[2] - point[2])
      
      # fix for identical 
      if(identical(point,p2) && x_len>y_len){ # this works for horizontal line
        xn_len <- line_len/2
        yn_len <- 0
      }else if(identical(point,p2) && x_len<y_len){ # this works for vertical line
        xn_len <- 0
        yn_len <- line_len/2
      }
      
      # Invert the x and y lengths and add/subtract from the refrence point
      # ref_points <- rbind(point,c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
      ref_points <- rbind(c(point[1] + yn_len,point[2] - xn_len),c(point[1] - yn_len,point[2] + xn_len))
      
      # use the reference points to return a line
      return(
        ref_points %>% 
          dplyr::as_tibble() %>% 
          dplyr::rename_with(tolower) %>% 
          sf::st_as_sf(coords = c("x", "y"), crs = my_crs, remove = F) %>%
          dplyr::summarise(n=dplyr::n()) %>%
          sf::st_cast("LINESTRING") %>%
          dplyr::ungroup() %>% 
          dplyr::select(-c(n))
      )
    }

function to pass a clump

generate data using clumping functions above, pass that data, return tree list with cut/keep flag based on the target clump group size

# treetops_clustered %>%
#   dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == prescription_target_data$min_clump_n_trees[3]) %>%
#   dplyr::arrange(desc(clump_n_trees)) %>%
#   # dplyr::slice(1) %>% 
#   # dplyr::pull(clump_id)
#   dplyr::filter(clump_id == 2) %>% 
#   dplyr::select(treeID, tidyselect::starts_with("clump_")) %>% 
#   dplyr::glimpse()
# levels(treetops_clustered$clump_n_trees_lb) %>% as.character() %>% as.numeric()
# c(
#   0 # always pad with 0
#   , levels(treetops_clustered$clump_n_trees_ub) %>% 
#     as.character() %>% 
#     as.numeric() %>%
#     unique() %>% 
#     purrr::discard(~ is.null(.x) || is.na(.x))
#   , Inf # always pad with Inf
# ) %>% 
# unique() %>% 
# sort() 
# function to cut a single clump
# as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb)))

# let's define a function to work with the target and the factors from st_clump_points
check_clump_breaks <- function(st_clump_points_ans) {
  # check cols
  cloud2trees:::check_df_cols_all_missing(
      st_clump_points_ans
      , col_names = c("clump_n_trees_grp","clump_n_trees_lb", "clump_n_trees_ub","clump_dist_m")
      , all_numeric = F, check_vals_missing = F
    )
  # clump distance used
  # get clump_dist_m used to create clumps in st_clump_points()
  clump_dist_m <- st_clump_points_ans %>% 
    dplyr::mutate(clump_dist_m = as.numeric(clump_dist_m)) %>% 
    dplyr::filter(!is.na(clump_dist_m)) %>% 
    dplyr::slice(1) %>% 
    dplyr::pull(clump_dist_m)
  if(length(clump_dist_m)!=1){stop("could not identify `clump_dist_m`")}
  
  # get lb and ub levels
  levels_clump_n_trees_lb <- levels(st_clump_points_ans$clump_n_trees_lb)
  levels_clump_n_trees_ub <- levels(st_clump_points_ans$clump_n_trees_ub)
  levels_clump_n_trees_grp <- levels(st_clump_points_ans$clump_n_trees_grp)
  # numeric
  num_clump_n_trees_lb <- levels_clump_n_trees_lb %>% as.character() %>% as.numeric()
  num_clump_n_trees_ub <- levels_clump_n_trees_ub %>% as.character() %>% as.numeric()
  
  ##### we need to id the cut_breaks used in st_clump_points() 
  cut_breaks <-
    c(
      0 # always pad with 0
      , num_clump_n_trees_ub %>% 
        unique() %>% 
        purrr::discard(~ is.null(.x) || is.na(.x))
      , Inf # always pad with Inf
    ) %>% 
    unique() %>% 
    sort() 
  if(identical(cut_breaks,c(0,Inf))){
    stop("could not identify cut breaks based on 'clump_n_trees_ub' factor levels")
  }
  
  return(list(
    levels_clump_n_trees_lb = levels_clump_n_trees_lb
    , levels_clump_n_trees_ub = levels_clump_n_trees_ub
    , levels_clump_n_trees_grp = levels_clump_n_trees_grp
    , num_clump_n_trees_lb = num_clump_n_trees_lb
    , num_clump_n_trees_ub = num_clump_n_trees_ub
    , cut_breaks = cut_breaks
    , clump_dist_m = clump_dist_m
  ))
  
}
# check_clump_breaks(treetops_clustered)

cut_clump_fn <- function(
  c # clump id from need_cut_trees data 
  , need_cut_trees # data with clumps already defined returned by st_clump_points()
  , tgt # based on: 
    # as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb))) from st_clump_points() 
    # -or- min_clump_n_trees from get_target_check_prescription()
){
  #### check target
  tgt <- as.numeric(as.character(tgt))
  if(!inherits(tgt,"numeric") || is.na(tgt)){
    stop("`tgt` does not meet numeric conditions")
  }
  # does this exist as level in the data?
  check_clump_breaks_ans <- check_clump_breaks(need_cut_trees)
  lvl_exists <- 
    check_clump_breaks_ans$num_clump_n_trees_lb %>% 
    purrr::has_element(tgt)
  if(!lvl_exists){
    stop("`tgt` does not exist in 'clump_n_trees_lb' factor of `need_cut_trees` data")
  }
  ##### we need to id the cut_breaks used in st_clump_points() to generate data in `need_cut_trees`
  these_cut_breaks <- check_clump_breaks_ans$cut_breaks
  # get distance used to create clumps in st_clump_points()
  dist_temp <- check_clump_breaks_ans$clump_dist_m
  
  # use only the current clump
  curr_dta <- need_cut_trees %>% 
    dplyr::filter(clump_id == c)
  
  if(nrow(curr_dta)<1){
    stop("cannot find data with the clump_id == `c` in `need_cut_trees` data")
  }
  
  # do we even need to cut?
  curr_clump_n_trees_lb <- curr_dta$clump_n_trees_lb %>% 
    as.character() %>% 
    as.numeric() %>% 
    unique()
  
  if(length(curr_clump_n_trees_lb)!=1){
    stop("multiple clumps identified in `need_cut_trees` data based on 'clump_n_trees_lb'")
  }
  # do we even need to cut?
  if(
    curr_clump_n_trees_lb<=tgt
    || nrow(curr_dta)<=tgt
  ){
    message("no need to cut as clump contains fewer trees than `tgt`")
    return(
      curr_dta %>% 
        sf::st_drop_geometry() %>% 
        dplyr::mutate(is_keep_tree=1) %>% 
        dplyr::select(treeID, is_keep_tree)
    )
  }
  
  # create clump polygon
  clumps <- get_clump_summary(curr_dta)
  
  # get the farthest line between points
  f_line_temp <- st_points_to_line(curr_dta, line_ext = dist_temp)
  
  # get perpendicular lines in dataset which we can iterate over to make cuts
  perp_line_sf_temp <-  
    # for every 1 m along line length, get a new perp line
    seq(
        from = 0
        , to = sf::st_length(f_line_temp) %>% as.numeric() %>% floor()
        , by = 1
      ) %>% 
    purrr::map(
      \(x)
      st_perp_line(
        interval=x, my_line = f_line_temp, proportion=F
      )
    ) %>% 
    dplyr::bind_rows() %>% 
    dplyr::mutate(line_n = dplyr::row_number())
  
  # find intersection of lines with the polygon and add length of intersection to perp line data
  perp_line_sf_temp <- perp_line_sf_temp %>% 
    dplyr::inner_join(
      # intersect and calc len
      perp_line_sf_temp %>% 
        sf::st_intersection(
          clumps %>% 
            dplyr::ungroup() %>% 
            dplyr::select(clump_id) %>% 
            dplyr::filter(clump_id == c)
        ) %>% 
        dplyr::mutate(len_m = sf::st_length(geometry) %>% as.numeric()) %>% 
        sf::st_drop_geometry()
      , by = "line_n"
    )
  
  # make cuts at the points where there is the least overlap with the clump polygon
    # list of potential line cut + tree combinations
    # aggregate the mean length of the intersecting cut lines to the tree level
    # prioritize trees for removal that have smallest length
    cut_tree_lines_temp <-  
      curr_dta %>% 
        sf::st_buffer(dist_temp/2) %>% 
        sf::st_join(perp_line_sf_temp) %>% 
        sf::st_drop_geometry() %>% 
        dplyr::group_by(treeID) %>% 
        dplyr::summarise(len_m = mean(len_m, na.rm = T)) %>% 
        dplyr::ungroup() %>% 
        dplyr::arrange(len_m, treeID) %>% 
        dplyr::mutate(n = dplyr::row_number())
  ###############################################################
  # while
  ###############################################################
  # cut until the desired clump size is achieved based on these cut lines
  while_temp <- 1
  i <- 1
  while(while_temp==1) {
    # cut trees
    cut_trees_temp <- cut_tree_lines_temp %>% 
      dplyr::slice(1:i) %>% 
      dplyr::distinct(treeID)
    
    # get the remaining trees not cut
    trees_remain_temp <- curr_dta %>% 
      dplyr::anti_join(cut_trees_temp, by = "treeID")
    
    # ensure that there are trees
    if(nrow(trees_remain_temp)==0){
      if(best_desired_grps_n_temp==0){
        # there are no more possible cuts :/
        # ... so we're going to add trees until the desired clump size is reached
        # start with biggest tree and keep adding trees until the desired clump size is reached
        if( any(stringr::str_equal(names(curr_dta), "dbh_cm")) ){
          start_tree <- curr_dta %>% 
            dplyr::arrange(desc(dbh_cm)) %>% 
            dplyr::slice(1)
        }else if( any(stringr::str_equal(names(curr_dta), "tree_height_m")) ){
          start_tree <- curr_dta %>% 
            dplyr::arrange(desc(tree_height_m)) %>% 
            dplyr::slice(1)
        }else{
          start_tree <- curr_dta %>% dplyr::slice(1)
        }
          
          increment_m <- 0.5
          k <- 1
          keep_trees <- dplyr::tibble(treeID = character(0))
          while(nrow(keep_trees)==0){
            # create a polygon to intersect with tree points and keep trees until number of trees met
            poly_keep <- start_tree %>% 
              sf::st_buffer( (dist_temp/2) + increment_m*k ) %>% 
              dplyr::mutate(dummy = 1) %>% 
              dplyr::select(dummy)
            # intersect and clump
            i_trees <- curr_dta %>% 
              sf::st_intersection(poly_keep) %>% 
              st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks)
            # huh?
            # i_trees %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp, clump_id)
            # keep only the desired clump
            keep_trees <- i_trees %>% 
              sf::st_drop_geometry() %>% 
              dplyr::filter(
                clump_id == 
                  i_trees %>% 
                    dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>% 
                    dplyr::pull(clump_id) %>% 
                    .[1] %>% 
                    dplyr::coalesce("nope")
              ) %>% 
              dplyr::select(treeID)
            # increment
            k <- k+1
          }
          # the remaining trees are cuts
          best_cuts <- curr_dta %>% 
            sf::st_drop_geometry() %>% 
            dplyr::anti_join(keep_trees, by = "treeID") %>% 
            dplyr::select(treeID)
      }else if(
        best_desired_grps_n_temp>0 
      ){ 
        ###############################
        # add trees back in until gets worse
        ###############################
        j <- 1
        while_add  <- 1
        while(while_add==1){
          # cut trees
          cut_trees_temp = best_cuts %>% 
            # add trees back in (i.e. remove them from the cut trees)
            dplyr::anti_join(
              cut_tree_lines_temp %>% 
                dplyr::slice(1:j) %>% 
                dplyr::distinct(treeID)
              , by = "treeID"
            )
          # get the remaining trees not cut
          trees_remain_temp <- curr_dta %>% 
            dplyr::anti_join(cut_trees_temp, by = "treeID")
          
          # count the groups remaining after cuts
          desired_grps_n_temp <- trees_remain_temp %>% 
            st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks) %>% 
            sf::st_drop_geometry() %>% 
            # do we have group sizes we want?
            dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>% 
            nrow()
          
          if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
            best_cuts <- cut_trees_temp
            best_desired_grps_n_temp <- desired_grps_n_temp
          }else{
            # stop it
            while_add <- 0  
          }
          # increment
          j <- j + 1
        } # while(while_add==1)
      } # best_desired_grps_n_temp>0 
      # done so stop the whole stop it
        while_temp <- 0
    }else{ # if(nrow(trees_remain_temp)==0)
      # count the groups remaining after cuts
      desired_grps_n_temp <- trees_remain_temp %>% 
        st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks) %>% 
        sf::st_drop_geometry() %>% 
        # do we have group sizes we want?
        dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>% 
        nrow()
      ### store best cut list
      if(i==1){
        best_cuts <- cut_trees_temp
        best_desired_grps_n_temp <- desired_grps_n_temp
      }else if(desired_grps_n_temp>best_desired_grps_n_temp){ # is this better than the best
        best_cuts <- cut_trees_temp
        best_desired_grps_n_temp <- desired_grps_n_temp
      }else if(
        desired_grps_n_temp==best_desired_grps_n_temp
        && i!=nrow(cut_tree_lines_temp)
      ){
        best_cuts <-  best_cuts
        best_desired_grps_n_temp <- best_desired_grps_n_temp
      }else if(
        best_desired_grps_n_temp>0 
        && desired_grps_n_temp<best_desired_grps_n_temp
      ){ # is this worse than the best which was successful?
        ###############################
        # add trees back in until gets worse
        ###############################
        j <- 1
        while_add <- 1
        while(while_add==1){
          # cut trees
          cut_trees_temp <- best_cuts %>% 
            # add trees back in (i.e. remove them from the cut trees)
            dplyr::anti_join(
              cut_tree_lines_temp %>% 
                dplyr::slice(1:j) %>% 
                dplyr::distinct(treeID)
              , by = "treeID"
            )
          # get the remaining trees not cut
          trees_remain_temp <- curr_dta %>% 
            dplyr::anti_join(cut_trees_temp, by = "treeID")
          
          # count the groups remaining after cuts
          desired_grps_n_temp <- trees_remain_temp %>% 
            st_clump_points(clump_dist_m = dist_temp, clump_breaks = these_cut_breaks) %>% 
            sf::st_drop_geometry() %>% 
            # do we have group sizes we want?
            dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == tgt) %>% 
            nrow()
          
          if(desired_grps_n_temp>=best_desired_grps_n_temp){ # is this better than the best
            best_cuts <- cut_trees_temp
            best_desired_grps_n_temp <- desired_grps_n_temp
          }else{
            # stop it
            while_add <- 0  
          }
          # increment
          j <- j + 1
        } # while(while_add==1)
        # done so stop the whole stop it
        while_temp <- 0
      }else if( i==nrow(cut_tree_lines_temp) ){ # is this the end?
        # stop it
        while_temp <- 0
      }
      ### increment
      i <- i+1
    } # else if(nrow(trees_remain_temp)==0)
  } # while(while_temp==1)
  
  # return it
  # return treelist with cut/keep
  # join to original data and pull
  d_temp <- curr_dta %>%
    sf::st_drop_geometry() %>%
    dplyr::mutate(is_keep_tree = as.numeric(NA)) %>% 
    dplyr::select(-c(is_keep_tree)) %>% 
    dplyr::left_join(
      best_cuts %>% dplyr::mutate(is_keep_tree = 0)
      , by = dplyr::join_by("treeID")
    ) %>%
    dplyr::mutate(is_keep_tree = dplyr::coalesce(is_keep_tree, 1)) %>%
    dplyr::select(treeID, is_keep_tree)
  # returns treeID and keep tree flag data frame
  return(d_temp)
}
# # example
# # pass it a clump id in the data
# check_clump_breaks_temp <- check_clump_breaks(treetops_clustered)
# # check_clump_breaks_temp
# xxx_temp <- treetops_clustered %>%
#   dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == check_clump_breaks_temp$num_clump_n_trees_lb[3] ) %>%
#   dplyr::arrange(desc(clump_n_trees)) %>%
#   dplyr::slice(1) %>%
#   # dplyr::select(tidyselect::starts_with("clump_"))
#   dplyr::pull(clump_id) %>%
#   cut_clump_fn(
#     need_cut_trees = treetops_clustered
#     , tgt = check_clump_breaks_temp$num_clump_n_trees_lb[2]
#   )
# dplyr::glimpse(xxx_temp)
# # now we have to look at the clumps
# zzz_temp <-
#   st_clump_points(
#     x = treetops_sf_with_dbh %>%
#       dplyr::inner_join(xxx_temp %>% dplyr::filter(is_keep_tree==1), by = "treeID")
#     , clump_dist_m = check_clump_breaks_temp$clump_dist_m
#     # , clump_breaks = these_cut_breaks
#     , clump_breaks = check_clump_breaks_temp$cut_breaks
#   )
# zzz_temp  %>%
#   dplyr::select(
#     treeID, tidyselect::starts_with("clump_")
#   ) %>%
#   dplyr::glimpse()
# zzz_temp  %>%
#   # rename the active geometry column to "geometry"
#   sf::st_set_geometry("geometry") %>%
#   dplyr::mutate(geometry = sf::st_buffer(geometry, dist = clump_dist_m/2)) %>%
#   ggplot2::ggplot() +
#   ggplot2::geom_sf(
#     mapping = ggplot2::aes(fill = clump_id)
#     , alpha = 0.8
#   ) +
#   ggplot2::theme_void()

there are some cases where the cut_clumps_fn may not return the desired tree clumps based on a single pass using the initial cut lines…iterate over the function until all clumps are in the desired clump size or smaller.

get_cut_clump <- function(
  c # clump id from need_cut_trees data 
  , need_cut_trees # data with clumps already defined returned by st_clump_points()
  , tgt # based on: 
    # as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb))) from st_clump_points() 
    # -or- min_clump_n_trees from get_target_check_prescription()
) {
  ### check
  check_clump_breaks_ans <- check_clump_breaks(need_cut_trees)
  # we need to id the cut_breaks used in st_clump_points() to generate data in `need_cut_trees`
  my_cut_breaks <- check_clump_breaks_ans$cut_breaks
  # get distance used to create clumps in st_clump_points()
  my_clump_dist_m <- check_clump_breaks_ans$clump_dist_m
  
  #### check target
  tgt <- as.numeric(as.character(tgt))
  if(!inherits(tgt,"numeric") || is.na(tgt)){
    stop("`tgt` does not meet numeric conditions")
  }
  # does this exist as level in the data?
  lvl_exists <- 
    check_clump_breaks_ans$num_clump_n_trees_lb %>% 
    purrr::has_element(tgt)
  if(!lvl_exists){
    stop("`tgt` does not exist in 'clump_n_trees_lb' factor of `need_cut_trees` data")
  }
  
  # what if we need to keep cutting because could not find a solution based on initial cut lines?
  cuts_first <- cut_clump_fn(
    c = c
    , need_cut_trees = need_cut_trees
    , tgt = tgt
  )
  # did we get the target?
  new_clumps <- 
    need_cut_trees %>% 
      dplyr::inner_join(
        cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
        , by = "treeID"
      ) %>%
      st_clump_points(clump_dist_m = my_clump_dist_m, clump_breaks = my_cut_breaks)
  # new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
  
  while(
    (
      new_clumps %>% dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) > tgt) %>% nrow()
    ) > 0
  ){
    # redo the cut clump
    # redo the cut clump
    cuts_again <- new_clumps %>% 
      dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) > tgt) %>%
      dplyr::pull(clump_id) %>% 
      unique() %>% 
      purrr::map(
        \(x)
        cut_clump_fn(c = x, need_cut_trees = new_clumps, tgt = tgt)
      ) %>% 
      dplyr::bind_rows() %>% 
      dplyr::rename(updt = is_keep_tree)
    
    # update the original cut data
    cuts_first <- cuts_first %>%
      dplyr::left_join(
        cuts_again
        , by = "treeID"
      ) %>% 
      dplyr::mutate(is_keep_tree = dplyr::coalesce(updt, is_keep_tree)) %>% 
      dplyr::select(-c(updt))
    
    # reset the new clumps
    new_clumps <-  
      need_cut_trees %>% 
        dplyr::inner_join(
          cuts_first %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
          , by = "treeID"
        ) %>%
        st_clump_points(clump_dist_m = my_clump_dist_m, clump_breaks = my_cut_breaks)
    # new_clumps %>% sf::st_drop_geometry() %>% dplyr::count(clump_n_trees_grp)
  } # end while
  return(cuts_first) 
}

function to pass a tree list

function to pass a whole tree list of sf point data

##############################################
# get cut keep tree flag
##############################################
get_keep_tree_flag <- function(
  x # x = sf point data
  , tgt # based on: 
    # as.numeric(as.character(levels(treetops_clustered$clump_n_trees_lb))) from st_clump_points() 
    # -or- min_clump_n_trees from get_target_check_prescription()
  , clump_dist_m = NA # size (radius) of the epsilon neighborhood = maximum distance between points to add to clump
  , clump_breaks = NA # where to break the clump size groups
  , clump_labels = NA # what to name the clump size groups
) {
  # MAKE CUTS INDEPENDENT OF CURRENT CLUMP GROUPING...
  # what if we pass a tree list with trees already cut? and thus, potentially different clump group sizes than is 
  # currently defined in clump_n_trees_grp?
  # 1) apply the clump grouping 
  # 2) only go through the cut alg for clumps that need cutting to the target
  # 3) append the trees in target clump or lower to the list at the end
  
  #### check target
  tgt <- as.numeric(as.character(tgt))
  if(!inherits(tgt,"numeric") || is.na(tgt)){
    stop("`tgt` does not meet numeric conditions")
  }
  
  #############################
  # 1) apply the clump grouping 
  #############################
  new_clump_trees <- st_clump_points(
    x = x
    , clump_dist_m = clump_dist_m
    , clump_breaks = clump_breaks
    , clump_labels = clump_labels
    # , ostory_ht_m = ostory_ht_m
    # , ostory_dbh_cm = ostory_dbh_cm
  )
  
  ### check
  check_clump_breaks_ans <- check_clump_breaks(new_clump_trees)
  # we need to id the cut_breaks used in st_clump_points() to generate data in `new_clump_trees`
  my_cut_breaks <- check_clump_breaks_ans$cut_breaks
  # get distance used to create clumps in st_clump_points()
  my_clump_dist_m <- check_clump_breaks_ans$clump_dist_m
  
  #### check target
  # does this exist as level in the data?
  lvl_exists <- 
    check_clump_breaks_ans$num_clump_n_trees_lb %>% 
    purrr::has_element(tgt)
  if(!lvl_exists){
    stop("`tgt` does not exist in 'clump_n_trees_lb' factor of `need_cut_trees` data")
  }
  
  # need cutting still
  need_cut_trees <- new_clump_trees %>% 
    dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) > tgt)
  
  #############################
  # 2) only go through the cut alg for clumps that need cutting to the target
  #############################
  # for each clump that still needs cutting, iterate over and return keep/cut flag
  compl_need_cut_trees <- need_cut_trees %>% 
    dplyr::pull(clump_id) %>% 
    unique() %>% 
    purrr::map(
      \(x)
      get_cut_clump(
        c = x
        , need_cut_trees = need_cut_trees
        , tgt = tgt
      )
    ) %>%
    dplyr::bind_rows()
  
  #############################
  # 4) append the trees in target clump or lower to the list at the end
  #############################
  compl_need_cut_trees <- compl_need_cut_trees %>% 
    dplyr::bind_rows(
      new_clump_trees %>% 
        sf::st_drop_geometry() %>% 
        dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) <= tgt) %>% 
        dplyr::mutate(is_keep_tree = 1) %>%
        dplyr::select(treeID, is_keep_tree)
    )
  
  # return
  return(
    # original data so that order is preserved
    x %>% 
      # add on trees that got cut with trees already good
      dplyr::inner_join(
        compl_need_cut_trees
        , by = "treeID"
      ) %>% 
      dplyr::pull(is_keep_tree)
  )
  
} # get_keep_tree_flag

Implement Cutting

This is currently very messy and needs to be compiled into an easily accessible function…;/

cut groups incrementally to achieve target tree clump group size proportions

must first have:

  • clustered trees from get_tree_clumps() named as treetops_clustered
  • clump polygons from get_clump_summary() named as clump_polys
  • target data from get_target_check_prescription() named as prescription_target_data
# must first have: 
# 
# * clustered trees from `get_tree_clumps()`
# * clump polygons from `get_clump_summary()`
# * target data from `get_target_check_prescription()`
dir_temp <- "../data/"
fnm_temp <- file.path(dir_temp,"prescription_trees.gpkg") # what should the compiled rgb be called?
if(!dir.exists(dir_temp)){dir.create(dir_temp, showWarnings = F)}
if(!file.exists(fnm_temp)){
  # check_clump_breaks(clump_polys) # will be the same
  check_clump_breaks_temp <- check_clump_breaks(treetops_clustered) # will be the same
  last_clump_n_trees_lb_temp <- min(check_clump_breaks_temp$num_clump_n_trees_lb)
  # what is the smallest group with a target?
  sm_grp_temp <-
    prescription_target_data %>% 
    dplyr::filter(pct_stand_n_trees_target>0) %>% 
    dplyr::pull(min_clump_n_trees) %>% 
    # as.numeric(as.character(clump_n_trees_lb))
    min()
  
  # start with the largest group in target data currently with trees
  grp_temp <-
    prescription_target_data %>% 
    dplyr::arrange(desc(min_clump_n_trees)) %>% 
    dplyr::filter(pct_stand_n_trees_current>0) %>% 
    dplyr::slice(1) %>% 
    dplyr::pull(min_clump_n_trees)
  
  # sm_grp_temp
  # grp_temp
  message(
    paste("started cutting", grp_temp, "at", Sys.time())
  )
  # identify clumps to leave untouched
  # clump_polys = get_clump_summary(ttops_temp)
  keep_clumps_temp <-
    clump_polys %>% 
      sf::st_drop_geometry() %>% 
      dplyr::mutate(
        min_clump_n_trees = as.numeric(as.character(clump_n_trees_lb))
      ) %>% 
      #keep only trees in current group
      dplyr::inner_join(
        prescription_target_data %>% 
          dplyr::filter(min_clump_n_trees == grp_temp) %>% 
          dplyr::select(
            mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
            , stand_n_clumps_target
          )
        , by = "min_clump_n_trees"
      ) %>% 
      # keep only clumps that meet criteria
      dplyr::filter(
        n_trees >= min_clump_n_trees
        & n_trees <= max_clump_n_trees
      ) %>% 
      # keep clumps closest the mean number in target
      dplyr::mutate(
        pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
      ) %>% 
      dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>% 
      dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>% 
      dplyr::select(clump_id)
  
  # dplyr::glimpse(keep_clumps_temp)
  
  # !!!!!!!!!!FIX: WHAT IF WE HAVE CLUMPS OF THIS SIZE BUT THE N_TREES>MAX_TREES AND NEED TO GET CLUMPS OF THIS SIZE?
  #   ... UPDATE CUTTING ALG TO CUT CLUMP DOWN TO MIN-MAX TREE RANGE FOR CLUMPS WITH INF UPPER LIMIT
  
  # start building tree list with keep/cut flag
  # ttops_temp = get_tree_clumps(
  #   my_suid = my_suid
  #   , clump_dist_m = clump_dist_m
  #   , ostory_dbh_cm = ostory_dbh_cm
  # )
  # build tree list 
  keep_trees <- 
    treetops_clustered %>% 
      dplyr::ungroup() %>% 
      sf::st_drop_geometry() %>% 
      dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>% 
      dplyr::select(treeID) %>% 
      dplyr::mutate(is_keep_tree = 1)
  
  # dplyr::glimpse(keep_trees)
  # did we keep all of the clumps?
  if( 
    nrow(keep_clumps_temp) == 
    ( clump_polys %>% dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == grp_temp) %>% nrow() )
  ){
    remaining_trees <- dplyr::tibble(
      treeID = character(0)
      , is_keep_tree = numeric(0)
    )
  }else{
    # determine keep/cut for the remaining trees in that group type 
    remaining_trees <- 
      treetops_clustered %>%
        dplyr::ungroup() %>% 
        # keep only trees in current grp
        dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) == grp_temp) %>% 
        # remove keep trees
        dplyr::anti_join(keep_trees, by = "treeID")
    
    # apply the cutting algorithm
    # get the flag
    remaining_trees$is_keep_tree <- get_keep_tree_flag(
        x = remaining_trees
        , tgt = prescription_target_data %>%  
          dplyr::ungroup() %>%
          dplyr::arrange(min_clump_n_trees) %>%
          dplyr::mutate(l = dplyr::lag(min_clump_n_trees)) %>%
          dplyr::filter(min_clump_n_trees == grp_temp) %>%
          dplyr::pull(l)
        , clump_dist_m = check_clump_breaks_temp$clump_dist_m
        , clump_breaks = check_clump_breaks_temp$cut_breaks
        , clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
      ) 
    
    # select relevant columns
    remaining_trees <- remaining_trees %>% 
      sf::st_drop_geometry() %>% 
      dplyr::ungroup() %>% 
      dplyr::select(treeID, is_keep_tree)
  }
  message(
    paste("done cutting", grp_temp, "at", Sys.time())
  )
  # remaining_trees %>% dplyr::count(is_keep_tree)
  
  ###############################################
  # now process to go on to the next groups
  ###############################################
  # get the next group
    grp_temp <-
      prescription_target_data %>%
      dplyr::filter(min_clump_n_trees<grp_temp) %>% 
      dplyr::arrange(desc(min_clump_n_trees)) %>% 
      dplyr::slice(1) %>% 
      dplyr::pull(min_clump_n_trees)
  
  while(grp_temp>=sm_grp_temp && grp_temp != last_clump_n_trees_lb_temp){
    message(
      paste("started cutting", grp_temp, "at", Sys.time())
    )
    # identify clumps to leave untouched
    # clump_polys = get_clump_summary(ttops_temp)
    keep_clumps_temp <- 
      clump_polys %>% 
        sf::st_drop_geometry() %>% 
        dplyr::mutate(
          min_clump_n_trees = as.numeric(as.character(clump_n_trees_lb))
        ) %>% 
        #keep only trees in current group
        dplyr::inner_join(
          prescription_target_data %>% 
            dplyr::filter(min_clump_n_trees == grp_temp) %>% 
            dplyr::select(
              mean_clump_n_trees, min_clump_n_trees, max_clump_n_trees
              , stand_n_clumps_target
            )
          , by = "min_clump_n_trees"
        ) %>% 
        # keep only clumps that meet criteria
        dplyr::filter(
          n_trees >= min_clump_n_trees
          & n_trees <= max_clump_n_trees
        ) %>% 
        # keep clumps closest the mean number in target group
        dplyr::mutate(
          pct_to_target = abs( (n_trees-mean_clump_n_trees)/mean_clump_n_trees )
        ) %>% 
        dplyr::arrange(pct_to_target, desc(qmd_cm), n_trees, clump_id) %>% 
        dplyr::filter(dplyr::row_number()<=stand_n_clumps_target) %>% 
        dplyr::select(clump_id)
    
    # add to tree list with keep/cut flag
    keep_trees <- keep_trees %>% 
      dplyr::bind_rows(
        treetops_clustered %>% 
          dplyr::ungroup() %>% 
          sf::st_drop_geometry() %>% 
          dplyr::inner_join(keep_clumps_temp, by = "clump_id") %>% 
          dplyr::select(treeID) %>% 
          dplyr::mutate(is_keep_tree = 1)
      )
  
    # check if the desired clump number was met and add the remaining trees from previous group if needed
    more_clumps_target <- 0
    if(
      nrow(keep_clumps_temp) <
      ( 
        prescription_target_data %>% 
          dplyr::filter(min_clump_n_trees == grp_temp) %>% 
          dplyr::pull(stand_n_clumps_target) 
      )
    ){
      # how many more clumps are needed?
      more_clumps_target <- 
        ( 
          prescription_target_data %>% 
            dplyr::filter(min_clump_n_trees == grp_temp) %>% 
            dplyr::pull(stand_n_clumps_target) 
        ) - nrow(keep_clumps_temp)
      
      # determine group size of remaining trees in the group prior that were not in a group selected and were not cut
      potential_trees <- 
        treetops_clustered %>% 
          dplyr::ungroup() %>% 
          dplyr::inner_join(
            remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
            , by = "treeID"
          ) %>% 
          st_clump_points(
            clump_dist_m = check_clump_breaks_temp$clump_dist_m
            , clump_breaks = check_clump_breaks_temp$cut_breaks
            , clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
          ) %>% 
          # ggplot() + geom_sf(aes(fill = clump_n_trees_grp)) + theme_void()
          # get the original clump id to prioritize new clump groups in the same area
          dplyr::inner_join(
            treetops_clustered %>% 
              sf::st_drop_geometry() %>% 
              dplyr::ungroup() %>% 
              dplyr::select(treeID, clump_id) %>% 
              dplyr::rename(orig_clump_id = clump_id)
            , by = "treeID"
          ) %>% 
          dplyr::group_by(orig_clump_id) %>% 
          dplyr::mutate(
            pct_desired_grp = 
              sum(ifelse( as.numeric(as.character(clump_n_trees_lb)) == grp_temp, 1, 0)) / dplyr::n()
          ) %>% 
          dplyr::ungroup() %>% 
          # keep only the current group
          dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == grp_temp)
      
      # pick trees from potential trees based on clump summary
      keep_trees <- keep_trees %>% 
        dplyr::bind_rows(
          potential_trees %>% 
            sf::st_drop_geometry() %>% 
            # filter trees based on clumps needed
            dplyr::inner_join(
              get_clump_summary(potential_trees) %>% 
                sf::st_drop_geometry() %>% 
                # join with original clump id metrics to prioritize 
                # keeping clumps in the same area and minimize cutting time
                dplyr::left_join(
                  potential_trees %>% 
                    sf::st_drop_geometry() %>% 
                    dplyr::group_by(clump_id, orig_clump_id) %>% 
                    dplyr::summarise(pct_desired_grp = max(pct_desired_grp))
                  , by = "clump_id"
                ) %>% 
                # keep the number of clumps needed
                dplyr::arrange(desc(pct_desired_grp), orig_clump_id, desc(qmd_cm), n_trees, clump_id) %>% 
                dplyr::filter(dplyr::row_number()<=more_clumps_target) %>%
                dplyr::select(clump_id)
              , by = "clump_id"
            ) %>% 
            dplyr::select(treeID) %>% 
            dplyr::mutate(is_keep_tree = 1)
        )
    } # end if don't have enough clumps
    ####################################################
    # determine keep/cut for the remaining trees in that group type and higher groups 
    ####################################################
    remaining_trees <- 
      treetops_clustered %>%
        dplyr::ungroup() %>% 
        # REMOVE TREES FROM PREVIOUS TREES REMAINING that got cut to make this group size
        dplyr::anti_join(
          remaining_trees %>% dplyr::filter(is_keep_tree == 0)
          , by = "treeID"
        ) %>%
        # keep only trees in current grp or prior grp
        dplyr::filter(as.numeric(as.character(clump_n_trees_lb)) >= grp_temp) %>% 
        # remove keep trees
        dplyr::anti_join(keep_trees, by = "treeID")
    
    if(nrow(remaining_trees) > 0){
      # apply the cutting algorithm
      # get the flag
      remaining_trees$is_keep_tree <- get_keep_tree_flag(
        x = remaining_trees
        , tgt = prescription_target_data %>%  
          dplyr::ungroup() %>%
          dplyr::arrange(min_clump_n_trees) %>%
          dplyr::mutate(l = dplyr::lag(min_clump_n_trees)) %>%
          dplyr::filter(min_clump_n_trees == grp_temp) %>%
          dplyr::pull(l)
        , clump_dist_m = check_clump_breaks_temp$clump_dist_m
        , clump_breaks = check_clump_breaks_temp$cut_breaks
        , clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
      ) 
      
      # select relevant columns
      remaining_trees <- remaining_trees %>% 
        sf::st_drop_geometry() %>% 
        dplyr::ungroup() %>% 
        dplyr::select(treeID, is_keep_tree)
    }
    
    # increment
      message(
        paste("done cutting", grp_temp, "at", Sys.time())
      )
    # get the next group
    grp_temp <-
      prescription_target_data %>%
      dplyr::filter(min_clump_n_trees<grp_temp) %>% 
      dplyr::arrange(desc(min_clump_n_trees)) %>% 
      dplyr::slice(1) %>% 
      dplyr::pull(min_clump_n_trees) %>% 
      dplyr::coalesce(last_clump_n_trees_lb_temp)
  }
  # now individuals
  if(sm_grp_temp == last_clump_n_trees_lb_temp){
    message(
      paste("started cutting", sm_grp_temp, "at", Sys.time())
    )
    potential_trees <- 
      # original data
      treetops_clustered %>% 
      dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == last_clump_n_trees_lb_temp) %>% 
      # remaining trees
      dplyr::bind_rows(
        treetops_clustered %>% 
          dplyr::ungroup() %>% 
          dplyr::inner_join(
            remaining_trees %>% dplyr::filter(is_keep_tree == 1) %>% dplyr::select(treeID)
            , by = "treeID"
          )
      ) %>% 
      # make sure that these are all individuals
      dplyr::mutate(is_orig = ifelse( as.numeric(as.character(clump_n_trees_lb))==last_clump_n_trees_lb_temp, 1, 0)) %>% 
      st_clump_points(
        clump_dist_m = check_clump_breaks_temp$clump_dist_m
        , clump_breaks = check_clump_breaks_temp$cut_breaks
        , clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
      ) %>% 
      dplyr::filter( as.numeric(as.character(clump_n_trees_lb)) == last_clump_n_trees_lb_temp)
    
      # pick trees from potential trees based on target
      keep_trees <- keep_trees %>% 
        dplyr::select(treeID) %>% 
        dplyr::bind_rows(
          potential_trees %>% 
            sf::st_drop_geometry() %>% 
            # keep the number of clumps needed
            dplyr::arrange(desc(is_orig), desc(dbh_cm), desc(tree_height_m)) %>% 
            dplyr::filter(
              dplyr::row_number() <= 
                (
                  prescription_target_data %>% 
                  dplyr::filter(min_clump_n_trees == last_clump_n_trees_lb_temp) %>% 
                  dplyr::pull(stand_n_clumps_target)
                )
            ) %>% 
            dplyr::select(treeID)
        ) %>% 
        dplyr::mutate(is_keep_tree = 1)
      
    message(
      paste("done cutting", sm_grp_temp, "at", Sys.time())
    )
  }
  # get the final prescription
  # treetops_clustered %>% dplyr::select(tidyselect::starts_with("clump_n_trees_")) %>% names()
  # treetops_clustered %>%
  #   dplyr::rename_with(
  #     .cols = tidyselect::starts_with("clump_")
  #     , .fn = ~ paste0("orig_", .x, recycle0 = TRUE)
  #   ) %>% 
  #   names()
  
  prescription_trees <- 
    treetops_clustered %>% 
    dplyr::ungroup() %>% 
    dplyr::left_join(
      keep_trees
      , by = "treeID"
    ) %>% 
    # tracking vars
    dplyr::mutate(
      # fill in keep tree flag
      is_keep_tree = dplyr::coalesce(is_keep_tree, 0)
    ) %>% 
    dplyr::rename_with(
      .cols = tidyselect::starts_with("clump_")
      , .fn = ~ paste0("orig_", .x, recycle0 = TRUE)
    )
  
  # attach the new clumping to the trees
  prescription_trees <- prescription_trees %>% 
    dplyr::left_join(
      # reclump
      prescription_trees %>% 
        dplyr::filter(is_keep_tree == 1) %>% 
        st_clump_points(
          clump_dist_m = check_clump_breaks_temp$clump_dist_m
          , clump_breaks = check_clump_breaks_temp$cut_breaks
          , clump_labels = check_clump_breaks_temp$levels_clump_n_trees_grp
        ) %>% 
        sf::st_drop_geometry() %>% 
        dplyr::select(treeID, tidyselect::starts_with("clump_")) %>% 
        dplyr::rename_with(
          .cols = tidyselect::starts_with("clump_")
          , .fn = ~ paste0("new_", .x, recycle0 = TRUE)
        )
      , by = "treeID"
    ) %>% 
    dplyr::mutate(
      new_clump_n_trees_grp = forcats::fct_na_value_to_level(new_clump_n_trees_grp, level = "Cut tree")
    )
  # save
  prescription_trees %>% sf::st_write(fnm_temp, append = F)
}else{
  prescription_trees <- sf::st_read(fnm_temp, quiet = T)
}

check the distribution of current vs prescription tree clump sizes

prescription_trees %>% 
  sf::st_drop_geometry() %>% 
  dplyr::count(
    orig_clump_n_trees_grp, new_clump_n_trees_grp
  ) %>% 
  tidyr::pivot_wider(names_from = orig_clump_n_trees_grp, values_from = n, values_fill = 0) %>% 
  dplyr::arrange(new_clump_n_trees_grp) %>% 
  dplyr::rename(`New` = new_clump_n_trees_grp) %>% 
  kableExtra::kbl() %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(
    c(" "=1,"Old"=length(unique(prescription_trees$orig_clump_n_trees_grp)))
  )
Old
New 10-15 trees 16-25 trees 2-4 trees 5-9 trees >25 trees Individual
10-15 trees 356 65 0 0 325 0
16-25 trees 0 35 0 0 139 0
2-4 trees 0 70 731 0 701 0
5-9 trees 0 38 0 607 514 0
>25 trees 0 0 0 0 228 0
Cut tree 0 184 0 0 3844 0
Individual 0 22 0 0 595 354

get the distance raster and openings vector data

treetops_clustered_new <-
  prescription_trees %>% 
  dplyr::filter(is_keep_tree == 1) %>% 
  dplyr::select(-tidyselect::starts_with("orig_")) %>% 
  # sf::st_drop_geometry() %>% dplyr::count(new_clump_n_trees_grp)
  dplyr::rename_with(
    .cols = tidyselect::starts_with("new_clump_")
    , .fn = ~ stringr::str_remove(.x, "^new_")
  )
# dplyr::glimpse(treetops_clustered_new)
clump_polys_new <- get_clump_summary(treetops_clustered_new)
# dplyr::glimpse(clump_polys_new)
# clump_polys_new %>% dplyr::select(clump_id) %>% plot()
get_clump_dist_rast_ans2 <- get_clump_dist_rast(get_clump_summary_ans = clump_polys_new, stand_sf = stand_sf)

check out the openings information

# what are the openings
get_clump_dist_rast_ans2$openings_vect$openining_area_m2 %>% summary()
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     0.04     0.04     1.92    71.53    17.96 11077.28
get_clump_dist_rast_ans2$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::summary()
## |---------|---------|---------|---------|=========================================                                          
##      layer       
##  Min.   : 0.000  
##  1st Qu.: 0.000  
##  Median : 0.714  
##  Mean   : 1.641  
##  3rd Qu.: 2.500  
##  Max.   :21.877  
##  NA's   :29973

plot the prescription

plt_fnl_12 <- 
  # ggplot2::ggplot() + 
  plt_stand +
    # distance
    ggplot2::geom_tile(
      data = get_clump_dist_rast_ans2$dist_rast %>% terra::aggregate(2, cores = 4) %>% terra::as.data.frame(xy = T) %>% dplyr::rename(f=3)
      , mapping = ggplot2::aes(x=x, y=y, fill = f)
    ) +
    harrypotter::scale_fill_hp(
      "mischief"
      , na.value = "transparent"
      , name = "distance to\nnearest tree (m)"
      , limits = c(0,22)
    ) +
    # openings
    ggplot2::geom_sf(
      data = get_clump_dist_rast_ans2$openings_vect
      , mapping = ggplot2::aes(color = openining_area_m2), fill = NA
    ) +
    ggplot2::scale_color_gradient(
      low = "gray77", high = "gray11"
      , labels = scales::comma_format(accuracy = 1)
      , name = latex2exp::TeX("opening\narea ($\\m^2$)")
      , limits = c(0,11100)
    ) +
    # clumps
    ggnewscale::new_scale_fill() +
    ggplot2::geom_sf(
      data = clump_polys_new
      , mapping = ggplot2::aes(fill = clump_n_trees_grp)
      , color = NA
    ) +
    harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) +
    # tree points
    ggplot2::geom_sf(data = treetops_clustered_new, color = "gray88", shape = ".") +
    ggplot2::theme_void() +
    ggplot2::theme(
      legend.position = "right"
      , legend.direction = "vertical"
      , plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
      , legend.title= ggplot2::element_text(size=8)
      , legend.text= ggplot2::element_text(size = 7)
    )
## |---------|---------|---------|---------|=========================================                                          

plot it

# plot
plt_fnl_12

# save it
ggplot2::ggsave("../data/06_clumps_dist_presciption.jpg", width = 7, height = 5)

highlight the openings

plt_open_12 <-
  plt_stand +
      # clumps
      ggplot2::geom_sf(
        data = clump_polys_new
        , mapping = ggplot2::aes(fill = clump_n_trees_grp)
        , color = NA
      ) +
      harrypotter::scale_fill_hp_d("always", name = "clump size", drop = F) + 
      # openings
      ggnewscale::new_scale_fill() +
      ggplot2::geom_sf(
        data = get_clump_dist_rast_ans2$openings_vect
        , mapping = ggplot2::aes(fill = openining_area_m2), color = NA
      ) +
      ggplot2::scale_fill_gradient(
        low = "gray77", high = "gray11"
        , labels = scales::comma_format(accuracy = 1)
        , name = latex2exp::TeX("opening\narea ($\\m^2$)")
        , limits = c(0,11100)
      ) +
      # tree points
      ggplot2::geom_sf(data = treetops_clustered_new, color = "gray88", shape = ".") +
      ggplot2::theme_void() +
      ggplot2::theme(
          legend.position = "right"
          , legend.direction = "vertical"
          , plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
          , legend.title= ggplot2::element_text(size=8)
          , legend.text= ggplot2::element_text(size = 7)
        )

# plot
plt_open_12

ggplot2::ggsave("../data/07_clumps_openings_presciption.jpg", width = 7, height = 5)

combine them with patchwork

plt_fnl_12 + (plt_open_12 + theme(legend.position = "none")) 

# &
#   ggplot2::theme(
#     plot.subtitle = ggplot2::element_text(size = 10, hjust = 0.5, face = "bold")
#     , legend.title=ggplot2::element_text(size=8)
#     , legend.text=ggplot2::element_text(size = 7)
#   )
ggplot2::ggsave("../data/08_clumps_combined_prescription.jpg", width = 10, height = 5.5)

combine old vs new in plot

(plt_fnl_1 + labs(subtitle = "Pre-Treatment")) +
  (plt_fnl_12  + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
  theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))

# save it
ggplot2::ggsave("../data/09_clumps_dist_compared.jpg", width = 10, height = 7)

combine old vs new openings highlight

(plt_open_1 + labs(subtitle = "Pre-Treatment")) +
  (plt_open_12  + labs(subtitle = "Post-Treatment") + theme(legend.position = "none")) &
  theme(plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold"), legend.title=element_text(size=8), legend.text=element_text(size = 7))

ggplot2::ggsave("../data/10_clumps_openings_compared.jpg", width = 10, height = 7)

table comparison

current vs target vs prescription

clump_n_trees_grp_summary2 <-
  get_clump_n_trees_grp_summary(
    get_tree_clumps_ans = treetops_clustered_new
    , get_clump_summary_ans = clump_polys_new
    , stand_sf = stand_sf
  ) 
# huh?
clump_n_trees_grp_summary2 %>% dplyr::glimpse()
## Rows: 6
## Columns: 33
## $ stand_area_ha               <dbl> 30.4578, 30.4578, 30.4578, 30.4578, 30.457…
## $ clump_n_trees_grp           <chr> "10-15 trees", "16-25 trees", "2-4 trees",…
## $ clump_n_trees_lb            <chr> "10", "16", "2", "5", "26", "1"
## $ clump_n_trees_ub            <chr> "15", "25", "4", "9", "Inf", "1"
## $ n_trees                     <int> 746, 174, 1502, 1159, 228, 971
## $ mean_dbh_cm                 <dbl> 27.14707, 28.91884, 27.77161, 27.56157, 27…
## $ mean_tree_height_m          <dbl> 14.71029, 16.23725, 15.14332, 14.96687, 15…
## $ loreys_height_m             <dbl> 15.96605, 17.23514, 16.21269, 16.05672, 16…
## $ basal_area_m2               <dbl> 44.93656, 11.72872, 93.86354, 71.45057, 14…
## $ clump_area_ha               <dbl> 1.631763, 0.376619, 3.562791, 2.597137, 0.…
## $ stand_n_clumps              <int> 63, 9, 513, 185, 7, 971
## $ trees_per_ha                <dbl> 24.492902, 5.712822, 49.314127, 38.052645,…
## $ basal_area_m2_per_ha        <dbl> 1.4753710, 0.3850810, 3.0817568, 2.3458871…
## $ qmd_cm                      <dbl> 27.69399, 29.29583, 28.20776, 28.01665, 28…
## $ stand_trees_per_ha          <dbl> 156.9384, 156.9384, 156.9384, 156.9384, 15…
## $ stand_basal_area_m2         <dbl> 303.453, 303.453, 303.453, 303.453, 303.45…
## $ stand_basal_area_m2_per_ha  <dbl> 9.963063, 9.963063, 9.963063, 9.963063, 9.…
## $ pct_stand_basal_area        <dbl> 0.14808408, 0.03865086, 0.30931822, 0.2354…
## $ pct_stand_n_trees           <dbl> 0.15606695, 0.03640167, 0.31422594, 0.2424…
## $ stand_qmd_cm                <dbl> 28.43065, 28.43065, 28.43065, 28.43065, 28…
## $ mean_dbh_in                 <dbl> 10.69594, 11.39402, 10.94202, 10.85926, 11…
## $ qmd_in                      <dbl> 10.91143, 11.54256, 11.11386, 11.03856, 11…
## $ stand_qmd_in                <dbl> 11.20168, 11.20168, 11.20168, 11.20168, 11…
## $ mean_tree_height_ft         <dbl> 48.24976, 53.25817, 49.67010, 49.09134, 50…
## $ loreys_height_ft            <dbl> 52.36863, 56.53127, 53.17761, 52.66603, 53…
## $ basal_area_ft2_per_ac       <dbl> 6.431142, 1.678568, 13.433378, 10.225722, …
## $ stand_basal_area_ft2_per_ac <dbl> 43.42899, 43.42899, 43.42899, 43.42899, 43…
## $ trees_per_ac                <dbl> 9.919625, 2.313693, 19.972221, 15.411321, …
## $ stand_trees_per_ac          <dbl> 63.56006, 63.56006, 63.56006, 63.56006, 63…
## $ stand_area_ac               <dbl> 75.26123, 75.26123, 75.26123, 75.26123, 75…
## $ clump_area_ac               <dbl> 4.0320853, 0.9306256, 8.8036559, 6.4175263…
## $ basal_area_ft2              <dbl> 483.6971, 126.2479, 1010.3472, 769.0939, 1…
## $ stand_basal_area_ft2        <dbl> 3266.368, 3266.368, 3266.368, 3266.368, 32…

table it to compare targets versus prescription

prescription_target_data %>% 
  dplyr::select(
    clump_n_trees_grp
    , tidyselect::starts_with("pct_stand_n_trees")
    , tidyselect::starts_with("trees_per_acre_")
  ) %>% 
  dplyr::left_join(
    clump_n_trees_grp_summary2 %>% 
      dplyr::select(clump_n_trees_grp, pct_stand_n_trees, trees_per_ac) %>% 
      dplyr::rename(
        pct_stand_n_trees_presc = pct_stand_n_trees
        , trees_per_acre_presc = trees_per_ac
      )
    , by = "clump_n_trees_grp"
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::where(is.numeric)
      , ~ dplyr::coalesce(.x,0)
    )
  ) %>% 
  dplyr::relocate(
    pct_stand_n_trees_presc
    , .after = pct_stand_n_trees_target  
  ) %>% 
  dplyr::relocate(
    trees_per_acre_presc
    , .after = trees_per_acre_target  
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_stand_n_trees")
      , ~ scales::percent(.x,accuracy = 1)
    )
  ) %>%
  kableExtra::kable(
    caption = paste0(
      "Current stand (ft<sup>2</sup> ac<sup>-1</sup>): "
      , clump_n_trees_grp_summary$stand_basal_area_ft2_per_ac[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand QMD (in): "
      , clump_n_trees_grp_summary$stand_qmd_in[1] %>% scales::number(accuracy = 0.1)
      , "<br>Current stand TPA: "
      , clump_n_trees_grp_summary$stand_trees_per_ac[1] %>% scales::number(accuracy = 1)
    )
    , escape = F
    , digits = 1
    , col.names = c(
      "", rep(c("current","target","prescription"),2)
    )
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::add_header_above(
    c(" "=1,"% Trees"=3, "TPA"=3)
  )
Current stand (ft2 ac-1): 77.1
Current stand QMD (in): 11.0
Current stand TPA: 117
% Trees
TPA
current target prescription current target prescription
Individual 4% 25% 20% 4.7 16.0 12.9
2-4 trees 8% 40% 31% 9.7 25.6 20.0
5-9 trees 7% 25% 24% 8.1 16.0 15.4
10-15 trees 4% 10% 16% 4.7 6.4 9.9
16-25 trees 5% 0% 4% 5.5 0.0 2.3
>25 trees 72% 0% 5% 84.4 0.0 3.0

check the pre- and post-treatment summary

# table it
clump_n_trees_grp_summary2 %>% 
  dplyr::mutate(presc = "post-treatment") %>% 
  dplyr::bind_rows(
    clump_n_trees_grp_summary %>% 
      dplyr::mutate(presc = "pre-treatment")
  ) %>% 
  dplyr::arrange(clump_n_trees_grp, desc(presc)) %>% 
  dplyr::select(
    clump_n_trees_grp, presc
    , n_trees
    , mean_dbh_in
    , qmd_in
    , mean_tree_height_ft
    , loreys_height_ft
    , trees_per_ac
    , basal_area_ft2_per_ac, pct_stand_basal_area, pct_stand_n_trees
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      .cols = c(pct_stand_basal_area, pct_stand_n_trees) 
      , .fns = ~ scales::percent(.x, accuracy = 1)
    )
  ) %>% 
  kableExtra::kbl(
    digits = 1
    , escape = F
    , caption = "Pre- and Post-treatment overstory tree clump summary"
    , col.names = c(
      ".", ""
      , "trees"
      , "mean<br>DBH (in)"
      , "QMD (in)"
      , "mean<br>Ht. (ft)"
      , "Loreys<br>Ht. (ft)"
      , "TPA"
      , "BA<br>ft<sup>2</sup> ac<sup>-1</sup>"
      , "%<br>stand BA"
      , "%<br>stand trees"
      )
  ) %>% 
  kableExtra::kable_styling() %>% 
  kableExtra::collapse_rows(columns = 1, valign = "top")
Pre- and Post-treatment overstory tree clump summary
. trees mean
DBH (in)
QMD (in) mean
Ht. (ft)
Loreys
Ht. (ft)
TPA BA
ft2 ac-1
%
stand BA
%
stand trees
10-15 trees pre-treatment 356 10.5 10.8 47.7 52.8 4.7 3.0 4% 4%
post-treatment 746 10.7 10.9 48.2 52.4 9.9 6.4 15% 16%
16-25 trees pre-treatment 414 10.5 10.7 46.9 50.8 5.5 3.5 4% 5%
post-treatment 174 11.4 11.5 53.3 56.5 2.3 1.7 4% 4%
2-4 trees pre-treatment 731 10.9 11.1 49.8 54.1 9.7 6.5 8% 8%
post-treatment 1502 10.9 11.1 49.7 53.2 20.0 13.4 31% 31%
5-9 trees pre-treatment 607 10.6 10.8 47.9 52.3 8.1 5.2 7% 7%
post-treatment 1159 10.9 11.0 49.1 52.7 15.4 10.2 24% 24%
>25 trees pre-treatment 6346 10.9 11.0 48.6 51.4 84.4 55.7 72% 72%
post-treatment 228 11.0 11.2 50.1 53.6 3.0 2.1 5% 5%
Individual pre-treatment 354 10.9 11.1 50.4 55.1 4.7 3.2 4% 4%
post-treatment 971 11.6 11.7 53.7 55.6 12.9 9.6 22% 20%